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Abstract 


The  objective  of  this  study  is  to  numerically  investigate  the  fundamental  roles  that  mo¬ 
mentum  and  vorticity  injections  play  in  suppressing  flow  separation  over  a  canonical  airfoil. 
Open-loop  control  of  separated,  incompressible  flow  over  a  NACA  0012  airfoil  at  Re  =  23, 000 
is  examined  through  large-eddy  simulations.  In  particular  two  conditions  are  considered:  (1) 
a  =  6°  -  shallow  stall  and  (2)  a  =  9°  -  deep  stall  for  applications  of  flow  control.  This 
study  does  not  attempt  to  replicate  a  specihc  actuator  but  aims  to  independently  introduce 
wall-normal  momentum  and  vorticity  flux  into  the  flow  through  model  boundary  conditions. 
We  hnd  that  the  modihcation  to  the  flow  held  can  be  captured  by  quantifying  both  the 
effects  of  wall-normal  momentum  (coefficient  of  momentum)  and  wall-normal  vorticity  (de¬ 
rived  coefficient  of  circulation),  by  considering  a  newly  dehned  total  input  parameter  (total 
coefficient).  The  inhuence  of  spanwise  spacing  is  also  examined  and  is  shown  that  the  total 
coefficient  accounts  for  spacing,  as  long  as  the  actuators  are  spaced  sufficiently  far  enough 
to  avoid  destructive  interference.  The  result  from  this  study  is  hoped  to  lead  to  a  general 
approach  for  quantifying  the  control  input  for  a  family  of  actuators. 

Moreover,  the  study  has  developed  advanced  analysis  techniques.  First,  the  capability 
to  perform  bi-global  stability  analysis  has  been  developed  and  validated,  which  can  serve 
as  a  basis  for  physics-based  active  how  control  guided  by  the  knowledge  of  hydrodynamic 
instabilities.  Second,  as  part  of  modeling  complex  unsteady  hows  in  general,  ehorts  in 
this  study  have  led  to  the  initial  development  of  a  novel  network-theoretic  approach  in 
quantifying  nonlinear  interactions  present  in  vortical  hows.  Dense  huid  how  graphs,  with 
vortices  as  nodes  and  induced  velocity  as  edge  weights  are  distilled  to  the  key  structures 
using  spectral  sparsihcation,  while  preserving  nonlinear  dynamics  and  invariants.  We  have 
also  been  able  to  quantify  two-dimensional  turbulence  as  a  weighted  scale-free  network  and 
evaluate  its  resilience  to  perturbations.  The  network-based  approach  to  analyze  interactions 
in  huid  hows  should  provide  a  refreshing  perspective  to  examine  a  wide  range  of  unsteady 
how  phenomena. 
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Chapter  1 
Introduction 


Separated  flow  is  a  major  cause  for  reduced  aerodynamic  performance  of  airfoils  at  high 
angles  of  attack,  which  leads  to  drag  increase  and  lift  decrease.  An  airfoil  at  a  high  incident 
requires  the  oncoming  flow  to  maneuver  the  leading  edge,  which  creates  an  adverse  pressure 
gradient  along  the  suction  surface,  which  ultimately  results  in  a  separated  boundary  layer. 
Such  behavior  of  the  boundary  layer  can  be  avoided  by  increasing  momentum  in  the  stream- 
wise  direction  of  the  flow  in  order  to  oppose  the  adverse  pressure  gradient.  Historically, 
flow  control  has  been  used  to  avert  separated  flow  and  the  accompanying  detrimental  effects 
(Chang,  1976).  Additional  momentum  can  be  introduced  into  the  boundary  layer  either 
directly  or  by  utilizing  free  stream  momentum  to  energize  the  boundary  layer  (Gad-el-Hak, 
2000a).  Directly  adding  momentum  to  the  boundary  layer,  quite  simply,  forces  the  bound¬ 
ary  layer  in  the  direction  of  the  free  stream.  Alternatively,  one  can  utilize  the  free  stream 
momentum  by  inducing  or  enhancing  the  mixing  between  the  free  stream  and  fluid  within 
the  boundary  layer  to  energize  the  flow  near  the  wall. 


1.1  Flow  Control 

As  a  means  to  prevent  or  delay  separation,  active  and  passive  flow  control  actuators  can  be 
utilized  to  introduce  perturbations  to  the  flow  field  (Gad-el-Hak,  2000a;  Cattafesta  et  ai, 
2008;  Joslin  &  Miller,  2009;  Cattafesta  &  Sheplak,  2011).  Active  flow  control  is  defined  by  the 
addition  of  external  energy  to  the  flow,  and  can  be  performed  with  a  large  assortment  of  flow 
control  actuators  (Cattafesta  &  Sheplak,  2011),  such  devices  include  steady  and  unsteady 
blowing/suction  (Lachmann,  1961;  Wu  et  al,  1998),  synthetic  jets  (Glezer  &  Amitay,  2002), 
plasma  actuators  (Corke  et  al,  2010),  and  vortex  generator  jets  (Compton  &  Johnston, 
1992;  Selby  et  al,  1992;  Zhang,  2003).  Passive  flow  control  devices  modify  the  flow  without 
external  energy  input  and  include,  leading  edge  modification  (Pedro  &  Kobayashi,  2008; 
Skillen  et  al,  2015),  vortex  generators  (VGs)  (Kerho  et  al,  1993;  Lin,  2002),  and  riblets 
(?Choi  et  al,  1993).  The  flow  control  actuators  listed  above  do  not  encompass  all  devices 
that  have  been  developed,  but  provides  an  idea  of  the  extent  of  the  variety  of  actuators  in 
use. 
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The  actuators  mentioned  above  have  been  implemented  in  different  scenarios  to  modify 
separated  flows.  Such  effort  has  relied  mostly  on  experiments  as  the  studies  have  often  been 
closely  related  to  actuator  development.  Seifert  &  Pack  (1999)  performed  experiments  with 
synthetic  jets  for  flow  over  NACA  0012  and  0015  airfoils  at  Reynolds  numbers  (1.5  x  10®  < 
Re  <  23.5  X  10®)  and  Mach  numbers  (0.2  <  M^o  <  0.65).  At  post-stall  angles  of  attack 
[a  >  8°),  the  reattached  flow  with  control  results  in  approximately  50%  increase  in  lift 
and  50%  decrease  in  drag.  More  recently,  Rathay  et  al.  (2014a, 6)  and  Seele  et  al.  (2012, 
2013)  increased  side  force  on  scaled  vertical  tail  stabilizer  using  synthetic  jet  actuators  and 
sweeping  jet  actuators.  The  former  studies  laid  the  foundation  for  a  successful  full-scale  flow 
control  on  a  vertical  tail  using  sweeping  jet  actuators  (Whalen  et  al.,  2015)  and  eventual 
implementation  on  the  Boeing  ecoDemonstrator  757.  In  another  study.  Little  et  al.  (2010) 
used  a  dielectric  barrier  discharge  plasma  actuator  to  modify  the  separation  due  to  the 
deflection  flap  of  a  high-lift  airfoil  at  Re  =  0(10®).  Effects  of  control,  which  include  delaying 
separation  as  well  as  lengthening  and  flattening  the  separated  region,  were  dependent  on 
the  waveform  of  the  control  input. In  the  same  spirit,  using  micro-vortex  generators,  Lin 
et  al.  (1994)  reduced  flap  separation  over  a  three-element  airfoil  in  high-lift  configuration. 
By  mitigating  separation,  the  lift-to-drag  ratio  is  doubled.  The  effectiveness  of  the  sub¬ 
boundary  layer  vortex  generators  is  highly  dependent  on  the  geometry,  spacing,  height, 
and  angle.  With  the  aforementioned  control  approaches,  flow  separation  is  mitigated  with 
different  actuators  but  often  leveraging  similar  control  mechanisms.  The  large  number  of 
studies  in  separation  control  have  given  emergence  to  different  types  of  non-dimensional 
parameters  (e.g.,  momentum  coefficient,  swirl  coefficient,  blowing  ratio)  used  to  report  the 
control  input  for  each  device. 

The  difficulty  in  analyzing  the  flow  through  experiments  is  caused  by  the  complex  three- 
dimensional  nature  of  separated  flow  under  the  influence  of  control  with  a  wide  range  in 
spatial  scales  being  present  in  the  flow  field.  For  the  same  computational  analysis  of  flow 
control  at  moderate  Reynolds  number  also  requires  substantial  resource  to  perform  a  siz¬ 
able  number  of  parameter  study  with  high  fidelity.  While  challenging  and  computationally 
expensive,  numerical  simulations  have  also  been  performed  to  study  flow  control  on  NACA 
airfoils.  In  addition  to  resolving  the  baseline  conditions  at  Reynolds  numbers  similar  to  ex¬ 
periments,  replicating  an  actuator  introduces  added  complexity  (Raju  et  al,  2009).  Earlier 
numerical  studies  have  examined  the  effectiveness  of  blowing/suction  (Wu  et  al.,  1998;  Huang 
et  al.,  2004)  and  vortex  generators  (Shan  et  al.,  2008).  More  recently,  high-fidelity,  three- 
dimensional  LES  of  complex  interactions  between  the  flow  over  the  airfoil  and  synthetic  jet 
actuation  have  been  performed  by  You  et  al.  (2008)  at  Re  =  896,  000  including  the  internal 
actuator  geometry.  These  types  of  complex  interactions  were  further  investigated  by  Abe 
et  al.  (2013)  at  Re  =  63,  000,  in  which  the  perturbations  with  different  spanwise  wavelengths 
were  considered  to  model  internal  effects  of  the  synthetic  jet.  Based  on  this  study,  Sato  et  al. 
(2015)  modeled  the  effects  of  plasma  actuators  and  performed  a  parametric  sweep  to  deter¬ 
mine  cases  for  separation  mitigation.  These  studies  highlighted  the  difficulty  in  replicating 
and  modeling  the  influence  of  specific  actuator  inputs. 

Regardless  of  the  selected  type  of  actuator,  we  can  view  the  flow  to  be  affected  by 
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means  of  mass,  momentum,  energy,  and  vorticity  injections.  Active  flow  control  actuators 
can  in  an  obvious  manner  add  the  aforementioned  forcing  inputs.  In  contrast,  passive  flow 
control  actuators  introduce  perturbations  in  response  to  the  flow  negotiating  the  actuator 
geometry.  Instead  of  simulating  flow  control  using  each  and  every  type  of  actuator,  we  aim 
to  understand  the  response  of  the  separated  flow  to  fundamental  control  inputs  that  these 
actuators  add.  The  present  investigation  is  of  the  influence  of  wall-normal  momentum  and 
vorticity  injections  for  separated  flow  over  a  NACA  0012  airfoil  at  9°  for  Re  =  23,  000.  The 
three-dimensional  LES  herein  sheds  light  on  the  separated  flow  physics  and  identify  the  key 
actuation  mechanism  for  reattaching  the  flow  at  high  angles  of  attack. 

1.2  Vortical  Interactions 

As  part  of  modeling  and  quantihcation  efforts  on  nonlinear  interactions  in  complex  fluid 
flows,  we  consider  novel  approaches.  In  particular,  we  extend  network  analysis  (Newman, 
2010)  and  graph  theory  (Bollobas,  1998)  to  canonical  examples  of  point- vortex  dynamics  and 
two-dimensional  turbulence,  as  an  initial  attempt  to  develop  these  network  based  techniques. 

In  the  held  of  huid  mechanics,  there  have  been  extensive  studies  performed  to  capture  the 
behavior  of  complex  huid  how.  Lagrangian  based  methods,  such  the  vortex  methods,  allow 
us  to  simulate  the  unsteady  huid  how  (Leonard,  1980;  Cottet  &  Koumoutsakos,  2000).  These 
methods  involve  modeling  of  huid  how  with  point  vortices,  vortex  sheets,  vortex  hlaments 
or  vortex  patches  (Sahman,  1992;  Cottet  &  Koumoutsakos,  2000).  Low-order  representation 
of  these  vortex  models  (Wang  &  Eldredge,  2013;  Heniati  et  al,  2014)  have  been  proposed 
in  recent  years.  The  evaluation  of  the  velocity  held  for  vortex  methods  rely  often  on  fast 
summation  methods  (Greengard  &  Rokhlin,  1987)  for  reduced  computational  time.  In  the 
present  work,  we  use  network-theoretic  approach  to  capture  the  behavior  of  vortical  elements 
involved  in  complex  huid  hows. 

As  the  computational  approaches  for  vortex  dynamics  are  being  developed,  there  are  also 
ongoing  ehorts  in  how  modeling.  Reduced-order  models  have  been  utilized  successfully  to 
describe  unsteady  incompressible  and  compressible  hows  (Rowley  et  al,  2004;  Noack  et  al, 
2005).  One  such  approach  is  to  utilize  Galerkin  projection  to  derive  reduced-order  models 
using  spatial  bases  such  as  the  Proper  Orthogonal  Decomposition  (POD)  modes  (Berkooz 
et  al,  1993;  Holmes  et  al.,  1996).  The  reduced-order  model  distills  the  inhnite-dimensional 
Navier-Stokes  equations  to  model  equations  with  state  variables  having  signihcantly  reduced 
dimensions.  In  the  present  work,  we  also  aim  to  capture  the  essential  physics  of  unsteady 
huid  how  but  not  by  reducing  the  dimension  of  the  state  variable.  Instead,  we  examine  the 
interaction  between  the  elements  of  the  state  variables  and  sparsify  the  interactions  utilizing 
a  network-theoretic  approach.  We  refer  to  the  dynamical  model  based  on  sparsifying  the 
interactions  as  sparsified  dynamics  in  this  paper. 

Network  describes  how  components  are  linked  to  one  another.  We  can  represent  the 
components  by  points  (nodes)  and  the  connections  by  lines  (edges)  through  a  mathematical 
abstraction.  The  structure  comprised  of  these  nodes  and  edges  is  called  a  graph,  which  has 
been  studied  in  detail  in  the  held  of  graph  theory  (Bollobas,  1998).  Network  analysis  is 
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concerned  with  the  study  of  graphs  as  well  as  the  interaction  and  evolution  of  the  variables 
of  interest  on  graphs  over  time  and  space  (Newman,  2010).  The  framework  developed  in 
network  analysis  and  graph  theory  can  provide  insights  into  how  the  structure  of  a  network 
can  influence  the  overall  dynamics  taking  place  on  the  network.  Network  analysis  and  graph 
theory  by  nature  are  very  fundamental  and  generic,  which  enable  them  to  impact  a  wide 
range  of  applications,  including  the  analysis  of  biological  and  social  networks,  study  of  traffic 
ffows,  and  design  of  robust  power  grids  (Newman,  2010).  Biologists  and  medical  scientists 
use  network  analysis  to  determine  how  electrical  signals  travel  inside  the  brain  and  how 
abnormality  in  the  brain  network  connections  can  affect  the  normal  functionality  (Duarte- 
Cavajalino  et  al,  2012;  Owen  et  ai,  2013). 

In  epidemiology,  researchers  model  the  outbreak  of  diseases  on  the  population  network. 
Locations  with  high  concentration  of  population,  such  as  airports,  stations,  schools,  and  hos¬ 
pitals,  can  be  represented  on  a  network  with  large  number  of  connections.  Identifying  such 
locations  is  especially  critical  when  containment  measures  are  designed  to  control  outbreaks 
of  HIV  (Morris,  1993),  SARS  (Lloyd-Smith  et  al,  2005),  and  Influenza  (Glass  et  al,  2006; 
Cauchemeza  et  al,  2011).  Each  of  these  diseases  have  different  dynamics  and  an  associ¬ 
ated  network  structure.  Once  the  high-risk  groups  and  areas  are  identihed,  network  analysis 
can  assist  in  designing  and  implementing  prevention  and  combat  strategies  in  the  most  swift 
manner  with  limited  resources  (Salathe  &  Jones,  2010;  Robinson  et  al,  2012).  Network  anal¬ 
ysis  can  also  reveal  how  a  group  of  people  are  socially  connected  to  one  another  and  examine 
how  subgroups  within  a  population  are  interlinked  in  a  complex  manner  (Porter  et  al,  2005). 
Moreover,  network  analysis  has  been  utilized  in  electrical  engineering  to  determine  the  volt¬ 
ages  and  currents  associated  with  electrical  circuits  via  graph  representations  (Chen,  2004). 
In  the  aforementioned  applications  of  network  analysis,  the  connections  between  people  or 
elements  are  highlighted. 

In  the  present  work,  we  consider  representing  vortices  with  nodes  and  the  interactions 
amongst  the  vortices  with  edges.  By  utilizing  the  network-theoretic  framework  to  study 
vortex  dynamics,  we  highlight  the  connections  (edges)  that  the  vortices  have  in  the  flow 
field.  Such  analysis  emphasizes  how  a  collection  of  vortices  influence  each  other  through 
a  causal  point  of  view  on  a  network  structure.  We  believe  that  the  present  study  can 
provide  an  alternative  tool  to  analyze  how  vortices  or  flow  structures  interact  in  the  flow 
held  and  support  the  development  of  interaction-based  models  to  capture  unsteady  vortex 
dynamics.  Furthermore,  we  consider  the  use  of  graph  sparsihcation  as  a  tool  for  sparsifying 
the  interaction  between  the  point  vortices.  These  models  keep  the  nodes  intact  and  reduce 
the  number  of  edges  maintaining  the  dimensionality  of  the  original  system  unlike  reduced- 
order  models.  The  removal  of  edges  can  drastically  reduce  computational  cost  to  model  the 
full  dynamical  behavior,  sharing  the  same  spirit  as  reduced-order  models. 

The  chaotic  motion  of  a  large  number  of  vortices  in  turbulent  hows  is  caused  by  the 
induced  velocities  of  the  vortices  themselves.  What  makes  turbulence  rich  and  complex  are 
the  vortical  interactions  in  the  how  held  that  take  place  over  a  wide  range  of  length  scales 
(Tennekes  &  Lumley,  1972;  Hinze,  1975;  Frisch,  1995;  Pope,  2000;  Davidson,  2004;  Lesieur, 
2008).  Thus,  complete  understanding  of  turbulence  has  remained  a  challenge  to  this  day 
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because  of  its  high-dimensionality,  multi-scale  interactions,  nonlinearity  and  the  resulting 
chaos.  Network  science  provides  an  alternative  view  of  complex  fluid  flows  in  terms  of  a 
network  of  vortex  interactions  (Nair  &  Taira,  2015),  and  this  perspective  illuminates  the 
underlying  structure  and  organization  of  turbulent  flows.  In  this  work,  we  show  that  two- 
dimensional  isotropic  turbulence  (Kraichnan  &  Montgomery,  1980;  McWilliams,  1984;  Benzi 
et  a/.,  1990;  Benzi  &  Colella,  1992;  Davidson,  2004;  Boffetta  &  Ecke,  2012)  has  a  scale-free 
network  structure  reminiscent  of  other  networks  found  in  nature  (Barabasi  &  Albert,  1999; 
Caldarelli,  2007).  While  most  of  the  attention  has  been  placed  on  unweighted  scale-free 
networks,  we  consider  the  use  of  weighted  scale-free  network  to  describe  the  variations  in  the 
strength  of  interactions  or  connectivities  (Barrat  et  ai,  2004).  Upon  identifying  the  network 
structure  of  turbulence,  physical  insights  can  be  obtained  as  to  which  vortical  interactions 
are  important  in  capturing  the  overall  physics  and  how  it  may  be  possible  to  control  the 
dynamics  of  turbulent  vortices  (Liu  et  al,  2011;  Farazmand  et  ai,  2011;  Brunton  &  Noack, 
2015). 

1.3  Structure  of  This  Report 

The  remainder  of  this  report  is  organized  as  follows.  We  present  the  computational  approach 
and  the  flow  control  setup  in  section  2.  Section  3  details  the  influence  of  control  on  the 
separated  flow,  with  focus  on  vortex  dynamics.  With  the  understanding  of  flow  control 
effects  on  the  flow,  we  quantify  control  inputs  in  section  4  and  are  related  to  the  change  in 
aerodynamic  forces.  The  flow  control  hndings  are  broadly  categorized  using  non-dimensional 
parameters  based  on  the  control  inputs.  We  perform  global  stability  analysis  in  section  5.  We 
highlight  network  analysis  methodology  in  section  6  and  quantify  the  vortical  interaction  in 
section  7.  We  extend  the  network  analysis  to  quantify  turbulent  flow  interactions  in  section 
8.  We  end  this  paper  with  concluding  remarks. 
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Chapter  2 

Flow  Control  Simulation  Approach 


2.1  Simulation  Setup 

We  perform  large-eddy  simulations  (LES)  of  flow  over  a  NACA  0012  airfoil  using  an  in¬ 
compressible  flow  solver,  Cliff  (CharLES  software  package,  Cascade  Technologies;  Ham  & 
laccarino  (2004);  Ham  et  al.  (2006)).  The  incompressible  Navier-Stokes  equations  are  dis¬ 
cretized  using  the  second-order  hnite-volume  and  time-integration  schemes.  Incorporating 
energy  conservation  properties  (Morinishi  et  ai,  1998),  the  solver  is  capable  of  handling 
structured  and  unstructured  grids.  In  the  present  computation,  we  dehne  the  Reynolds 
number  as  Re  =  Uooc/i'  =  23,000,  which  is  based  on  the  free  stream  velocity,  f/oo,  the 
chord  length  of  the  airfoil,  c,  and  the  free  stream  viscosity,  v.  To  predict  the  flow  held  at 
this  Reynolds  number,  eddy  viscosity  is  introduced  with  the  Vreman  subgrid-scale  model 
(Vreman,  2004).  We  investigate  the  inhuence  of  control  of  fully  separated  how  at  a  =  9°  for 
the  purpose  of  increasing  lift  and  reducing  drag. 

The  overall  size  of  the  computational  domain  is  {x / c^y / c,  z / c)  G  [—20,25]  x  [—20,20]  x 
[—0.1,  0.1],  as  illustrated  in  Fig.  2.1.  The  spatial  directions  x,  y,  and  2;  refer  to  the  streamwise, 
wall-normal,  and  spanwise  directions,  respectively.  At  the  inlet,  laminar,  uniform  how  of 
u/Uoo  =  (1,0,0)  is  prescribed.  Stress-free  boundary  conditions  are  applied  at  the  far-held 
(top  and  bottom)  boundaries.  A  convective  outhow  condition,  ^  =0,  is  prescribed 

at  the  outlet,  with  the  convective  velocity  {Uc)  set  to  the  time-averaged  normal  velocity,  to 
allow  wake  structures  to  leave  the  domain  without  disturbing  the  near-held  solution.  In  the 
spanwise  direction,  how  is  taken  to  be  periodic. 

The  present  study  utilizes  a  hybrid  structured/unstructured  spatial  discretization.  A 
structured  grid  is  used  to  achieve  adequate  resolution  in  the  near  held  of  the  airfoil  and  an 
unstructured  far-held  grid  utilized  to  reduce  the  number  of  points  needed  in  this  region  of  the 
computational  domain.  The  mesh  is  planar  (two-dimensional)  and  extruded  in  the  spanwise 
direction.  The  computational  domain  in  the  vicinity  of  the  actuators  is  also  rehned  in  order 
to  resolve  the  hne-scale  how  structures  produced  by  interaction  of  the  incoming  how  and  the 
actuator  inputs.  Based  on  the  domain  size  and  necessary  spatial  resolution,  the  resulting 
computational  domain  is  comprised  of  approximately  40  x  10®  grid  points. 
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Figure  2.1:  Computational  domain  used  for  the  LES,  with  the  background  flow  directed 
from  left  to  right.  A  top  down  view  to  the  right  shows  the  airfoil  with  two  actuators  over 
the  spanwise  extent. 


To  perform  wall-resolved  LES  in  the  present  investigation,  sufficient  resolution  is  needed 
near  the  airfoil  surface.  The  largest  non-dimensional  wall  spacing  along  the  suction  surface  of 
the  airfoil  is  x~^  6,  i/’*'  ~  0.5,  and  z~^  ~  12  (x~^  =  UrXjv,  where  Ur  is  the  friction  velocity). 

No  slip,  no  penetration  boundary  conditions  are  applied  on  the  airfoil  surface.  The  surface 
mesh  of  the  airfoil  is  rehned  such  that  30  points  span  across  each  actuator,  which  we  describe 
in  detail  below  in  Section  2.2.  For  the  controlled  cases,  an  actuator  velocity  prohle,  to  be 
discussed  later,  is  specihed  at  the  actuator  locations  on  the  surface  of  the  airfoil.  As  shown 
in  Fig.  2.1,  the  spanwise  actuator  spacing  is  set  to  l^/c  =  0.1  for  the  majority  of  the  present 
study.  Later  in  this  paper,  we  consider  using  1  to  4  actuators  {Iz/c  =  0.2  to  0.05)  over  the 
span  of  0.2  to  examine  the  influence  of  spanwise  spacing  on  control  effectiveness. 

Throughout  this  paper,  the  lift,  drag,  and  pressure  coefficients  are  dehned  as 


Cr  = 


Ft 


IpooUIA^ 


Cd  = 


2  Poo 


UlA^ 


and  Cp  = 


P- Poo 

In  m  ■ 

2^00'^  oo 


(2,1) 


where  A  =  l^c  is  the  planform  area  of  the  airfoil,  and  poo  is  the  free  stream  density.  The 
forces,  Fl  and  Fd,  represent  lift  and  drag,  respectively.  The  coefficient  of  pressure  subtracts 
off  the  free-stream  pressure  value  Poo  and  is  normalized  by  the  dynamic  pressure.  The  forces 
produced  by  control  are  included  in  the  reported  lift  and  drag  values. 

To  ensure  the  validity  of  our  computational  approach,  we  compare  our  computational 
results  for  the  baseline  flow  with  those  reported  in  Kojima  et  al.  (2013).  Upon  examining 
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the  flow  fleld,  airfoil  pressure  distribution,  and  forces,  we  observe  good  agreement  as  they 
have  been  reported  in  detail  in  our  prior  publications  (Munday  &  Taira,  2014,  2015).  More¬ 
over,  the  computational  requirements  for  the  computational  domain  size,  boundary  condition 
treatment,  and  grid  resolution  were  examined  with  care.  With  the  separated  flow  validated, 
we  now  consider  the  applications  of  flow  control  to  mitigate  separation. 

2.2  Actuation  Setup 

We  introduce  steady  circular  jets  with  swirl  on  the  suction  surface  of  the  airfoil  specified 
by  velocity  boundary  conditions.  A  canonical  setup  for  the  actuators  is  shown  in  Fig.  2.1 
(right).  The  actuator  jets  are  specified  by  circular  regions  with  a  radius  of  tq/c  =  0.01  and 
placed  at  10%  chord  location.  This  is  the  typical  actuator  location  used  in  the  majority  of 
the  cases,  but  we  also  report  the  influence  of  the  spanwise  spacing  and  chordwise  placement 
(herein  the  chordwise  direction  is  defined  as  ^).  In  particular,  we  consider  1  to  4  actuators 
across  the  spanwise  extent.  In  order  to  assess  the  influence  of  the  injection  of  wall-normal 
momentum  and  wall-normal  vorticity  independently,  wall-normal  (m„)  and  azimuthal  velocity 
(ug)  profiles  are  specified.  At  the  actuator  locations,  the  time-invariant  velocity  profiles  are 


(2.2) 


which  are  plotted  in  Fig.  2.2.  For  simulations  with  multiple  actuators  over  the  spanwise 
extent  of  the  computational  domain,  we  consider  arrangements  where  swirl  is  introduced  in 
a  co-rotating  or  counter-rotating  manner.  Throughout  the  paper,  we  denote  co-rotating  cases 
with  COR  and  counter-rotating  ones  with  CTR,  respectively.  The  magnitude  of  wall-normal 
velocity  is  chosen  such  that  the  coefficient  of  momentum  is  on  the  same  order  of  magnitude 
as  in  previous  successful  studies  (Greenblatt  et  ai,  2015).  The  azimuthal  velocity  is  selected 
such  that  the  maximum  velocity  is  on  the  same  order  of  magnitude  as  the  free  stream,  which 
is  the  upper  bound  for  vortex  generators. 

To  examine  the  fundamental  influence  of  the  above  actuation  inputs  on  separated  flow 
over  the  airfoil,  we  consider  a  large  number  of  flow  control  cases.  We  considered  control  at 
a  =  6°  and  a  =  9°  and  a  summary  of  all  the  controlled  cases  are  tabulated  in  Appendix 
A.l  and  A. 2,  respectively.  In  the  tables,  all  of  the  control  input  values  are  listed  along  with 
the  resulting  force  coefficients.  Also  listed  in  the  appendix  is  the  total  coefficient  which  will 
be  defined  and  discussed  in  detail  in  Section  4.  In  what  follows  our  discussion  pertains  to 
a  =  9°.  In  the  present  study,  we  keep  our  actuation  inputs  to  be  steady  to  limit  the  numbers 
of  parameters  and  LES  computations  to  be  tractable. 
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Figure  2.2:  Prescribed  velocity  profiles  for  flow  control  inputs. 
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Chapter  3 
Flow  Physics 


In  this  section,  we  examine  the  flow  fleld  for  both  nncontrolled  and  controlled  cases  in  detail. 
At  an  angle  of  attack  of  a  =  9°  and  Re  =  23,  000,  flow  separates  near  the  leading  edge  and 
does  not  reattach  which  is  consistent  with  results  of  similar  studies  (Chang,  1976;  Kotapati 
et  al,  2010;  Kojinia  et  al,  2013).  First,  we  discuss  how  the  separated  flow  will  be  evaluated 


3.1  Flow  Visualization 


To  highlight  the  effect  of  control  input  on  the  flow  fleld,  the  time-averaged  zero-streamwise- 
velocity  [ux  =  0,  herein)  iso-surface  is  visualized.  By  definition,  separated  flow  is  accom¬ 
panied  by  reverse  flow.  The  Ux  =  0  iso-surface  encompasses  the  reverse  flow  region  and 
is  used  here  to  capture  separation  and  reattachment.  While  the  u  =  0  surface  provides  a 
general  sense  of  the  size  of  the  separated  region,  it  does  not  encompass  and  should  not  be 
considered  as  the  entirety  of  a  separated  boundary  layer.  We  also  superpose  the  spanwise 
Reynolds  stress  distribution  on  the  Ux  =  0  iso-surface  to  provide  insights  into  turbulent 
mixing  of  near-wall  and  free  stream  momentum.  The  spanwise  Reynolds  stress  is  defined 
as  Txy  =  pu'^u'y,  which  is  based  on  the  fluctuations  in  the  x  and  ^/-directions,  u'^  and  u'y, 
respectively.  Visualizing  the  Reynolds  stress  identifies  locations  in  which  turbulent  mixing 
is  enhanced  to  aid  in  reattachment  with  control  input.  Figure  3.1,  illustrates  the  Ux  =  ^ 
iso-surface  for  the  baseline  case.  The  images  in  Fig.  3.1  are  of  the  time-averaged  flow  fields, 
and  the  red  line  highlights  the  locations  of  the  Ua,  =  0  iso-surface  amongst  the  time-averaged 
streamlines.  For  the  baseline  flow,  the  time-averaged  results  are  homogeneous  throughout 
the  spanwise  direction.  Identifying  locations  of  reversed  flow  throughout  the  span  will  help 
in  understanding  the  effects  of  control  input  later. 

For  incompressible  flow,  vorticity  can  only  be  generated  along  the  wall  surface  or  injected 
by  the  control  inputs.  Vorticity  is  introduced  to  the  flow  through  a  wall-normal  diffusive 
flux  at  the  surface 


z/(Vn;)o  •  n  =  —n  x 


du  1  „ 

— +  -Vp 
dt  p 


(3,1) 


J  0 


which  is  caused  by  the  acceleration  of  the  wall  and  the  local  pressure  gradient  in  the  tangential 
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Figure  3.1:  Time- averaged  streamlines  (black)  and  the  zero-streamwise-velocity  =  0) 
curve  (red)  on  the  left.  The  zero-streamwise-velocity  iso-surface  colored  with  Reynolds  stress 
is  to  the  right. 


direction  (Hornung,  1989;  Wu  et  al.,  2006).  Here,  the  subscript  zero  denotes  the  surface  value 
and  n  represents  the  wall-normal  unit  vector.  Since  the  wing  is  assumed  stationary  in  the 
present  study,  only  the  pressure  gradient  attributes  to  the  above  flux.  This  wall-normal 
diffusive  flux  is  referred  to  as  the  source  of  vorticity  (Hornung,  1989)  and  is  denoted  by  cr  ■  h 

(T  ■n=  =  v{Vu:)o-h,  (3.2) 

The  subscripts  of  CTj  refer  to  the  directions  of  the  wall-normal  flux  components.  Throughout 
the  study,  time-averaged  cxj  are  shown  to  closely  examine  the  local  influx  and  efflux  of 
vorticity  in  the  separated  region  of  the  flow. 

The  spanwise  vorticity  flux  for  the  baseline  case  is  shown  in  Fig.  3.2.  Negative  spanwise 
vorticity  flux  is  observed  upstream  of  the  separation  location  as  the  boundary  layer  negotiates 
the  leading  edge.  Flow  maneuvering  the  curvature  of  the  leading  edge  causes  a  pressure 
gradient  in  the  streamwise  direction,  which  is  correlated  to  regions  with  increased  spanwise 
vorticity  flux.  Additionally,  near  the  reattachment  location,  an  increase  in  spanwise  vorticity 
flux  in  the  same  direction  (seen  in  red)  is  observed.  The  spanwise  vorticity  flux  at  the 
reattachment  location  is  introduced  due  to  the  pressure  gradient,  as  shown  in  Fig.  3.2.  Time- 
averaged  quantities  are  used,  therefore  the  streamwise  and  wall-normal  vorticity  flux  are  zero 
for  the  baseline  cases.  Through  this  analysis,  emphasis  will  be  placed  on  the  separation  and 
reattachment  location  (if  present),  as  well  as  near  the  actuators  to  observe  the  effects  control 
has  on  the  vorticity  flux. 

3.2  Separation  Control 

Baseline  Flow 

Now,  let  us  consider  the  separated  flow  over  the  NACA  0012  airfoil.  Flow  over  the  airfoil  at 
this  higher  angle  of  attack  separates  from  the  leading  edge  and  does  not  reattach  over  the 
chord  as  seen  in  Fig.  3.3  (top).  A  shear  layer  develops  at  the  leading  edge  and  generates 
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Figure  3.2:  Spanwise  vorticity  flux  on  the  suction  surface  of  the  airfoil  for  the  baseline  flow 
case.  Surface  pressure  prohles  for  the  corresponding  cases  are  shown  on  the  right. 


Case 

"^n^max  /  U QiQ 

c. 

Cr 

Rot.  dir. 

Cd 

Cl 

Baseline 

0 

0 

0 

0 

- 

0.118 

0.584 

9A 

1.26 

0 

0.44% 

0 

- 

0.108 

0.403 

9B 

1.26 

0.95 

0.44% 

2.29% 

CTR 

0.094 

0.483 

9C 

1.26 

1.26 

0.44% 

2.60% 

COR 

0.065 

0.671 

Table  3.1:  Representative  flow  control  cases  for  a  =  9°  chosen  for  flow  visualization  in 
Fig.  3.3. 

large  spanwise  vortices.  The  breakdown  of  these  vortices  leads  to  turbulent  flow  past  the 
mid-chord  of  the  airfoil.  The  Ux  =  ^  iso-surface  covers  the  entire  suction  surface  of  the 
airfoil,  even  extending  downstream  into  the  wake.  The  resulting  large  region  of  reverse  flow 
is  detrimental  to  the  aerodynamic  forces  on  the  airfoil.  At  the  trailing  edge  bluff  body 
shedding  is  observed,  which  is  noted  by  the  large-scale,  opposite  sign  vortices.  Spanwise 
Reynolds  stress  on  the  Wa,  =  0  iso-surface  shows  the  development  of  a  large  turbulent  wake 
with  the  breakdown  of  the  spanwise  vortex.  Through  the  investigation  of  the  vorticity  flux  on 
the  surface  of  the  airfoil.  A  high  level  of  spanwise  vorticity  flux  is  observed  near  the  leading 
edge  and  is  correlated  to  a  strong  pressure  gradient  caused  by  the  curvature  of  the  airfoil. 
Since  there  is  no  reattachment,  no  other  large  influxes  of  vorticity  are  observed  downstream 
of  the  separation  location. 

In  an  attempt  to  reattach  the  massively  separated  flow  over  the  airfoil  at  a  =  9°,  actu¬ 
ation  with  wall- normal  velocity  of  Wn.max/f^oo  =  0  to  1.83  and  azimuthal  velocity  values  of 
ug^max/Uoo  =  0  to  2.52  is  considered.  The  actuator  location  is  based  on  the  separation  point 
for  the  a  =  6°  case  {^/c  =  0.1).  For  the  angle  of  attack  of  a  =  9°,  the  actuator  is  located 
downstream  of  the  separation  point  (we  later  examine  the  effect  of  moving  the  actuator  to 
the  separation  point  at  ^/c  =  0.05).  Three  representative  control  cases  for  a  =  9°  are  listed 
in  Table  3.1,  and  we  note  that  the  aerodynamic  forces  can  be  altered  signihcantly,  up  to  36% 
reduction  in  drag  and  31%  increase  in  lift  with  appropriate  control  inputs. 

With  only  pure  blowing  (case  9A,  Un,max/Uoo  =  1.26)  flow  control  negatively  impacts 
lift.  In  Fig.  3.3  we  observe  that  the  separated  flow  is  not  modihed  much  by  this  control 
conhguration.  Comparing  case  9A  to  the  baseline  flow,  the  size  of  the  separated  region 
increases  in  the  wake.  Although  a  dehcit  in  the  reverse  flow  is  created  in  the  vicinity  of  the 
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actuator  for  case  9A,  flow  eventually  separates  throughout  the  span.  We  found  solely  adding 
momentum  with  a  wall-normal  velocity  of  Mn,max/f^oo  =  1-26  is  not  effective  in  overcoming 
the  adverse  pressure  gradient  throughout  the  span/chord.  The  most  notable  difference  in 
the  instantaneous  flow  profiles  of  the  pure  blowing  case  is  the  spanwise  vortex  is  no  longer 
present.  The  actuator  jets  show  the  Kelvin-Helmholtz  instability  which  interacts  with  the 
separated  shear  layer.  This  results  in  the  spanwise  vortex  break  down  close  to  the  leading 
edge.  The  Reynolds  stress  on  the  iso-surface  is  of  lower  level  than  that  observed  for  the 
baseline  case.  The  reduced  amount  of  turbulent  mixing  from  the  actuator  perturbation  does 
not  result  in  reattaching  flow.  We  observe  secondary  surface  vorticity  flux  with  case  9A 
due  to  obstruction  and  resulting  pressure  gradients  created  by  the  actuator.  Due  to  the  flow 
separating  upstream  of  the  actuator,  flow  is  not  forced  between  the  actuators  as  drastically  as 
the  lower  angle  of  attack,  which  increased  momentum  between  the  actuators  for  the  previous 
angle  of  attack. 

Next,  let  us  introduce  swirling  motion  to  the  control  jets.  As  shown  in  Fig.  3.3,  the 
combination  of  wall-normal  momentum  and  vorticity  injections  (cases  9B  and  9C)  can  greatly 
modiiy  the  time-averaged  flow  field.  Cases  9B  and  9C  use  azimuthal  actuator  velocity  profiles 
with  Un,max/woo  =  0.95  and  1.26,  respectively.  For  these  two  cases,  the  addition  of  wall- 
normal  vorticity  and  wall-normal  momentum  decrease  the  size  of  the  reverse  flow  region. 
For  case  9B,  flow  still  separates  near  the  leading  edge  of  the  airfoil  and  reattaches  in  the 
vicinity  of  the  actuator,  forcing  the  oncoming  flow  to  be  redirected  closer  to  the  airfoil. 
Overall,  the  resulting  size  of  the  separated  flow  region  is  noticeably  smaller.  The  change 
in  the  time-averaged  flow  held  is  rehected  in  the  forces;  drag  is  further  decreased  and  lift 
increased  compared  to  pure  blowing  with  identical  wall-normal  velocity  input.  Once  again, 
the  large  structures  emanating  from  the  actuators  seen  in  the  pure  blowing  case  are  broken 
into  smaller  structures  when  rotation  is  added  to  actuation. 

Fully  reattaching  the  how  with  the  same  wall-normal  momentum  requires  additional  wall- 
normal  vorticity  input.  The  control  input  with  Mn,max/woo  =  1-26  and  Me,max/woo  =  1-26  by 
case  90  results  in  a  fully  attached  boundary  layer.  This  fully  attached  case  translates  to 
signihcant  improvements  in  terms  of  aerodynamic  forces.  In  case  90  we  observe  in  Fig.  3.3 
(bottom)  the  reverse  how  region  is  diminished  downstream  of  the  actuators.  Although  how 
separates  upstream  of  the  actuators,  but  the  momentum  added  to  the  boundary  layer  by  the 
actuators  allows  for  the  how  to  overcome  the  adverse  pressure  gradient.  The  wall-normal 
momentum  penetrates  the  boundary  layer  and  transports  injected  vorticity  into  the  shear 
layer,  mixing  the  free  stream  momentum  and  boundary  layer.  Minimal  Reynolds  stress 
is  observed  on  the  Ua,  =  0  iso-surface,  which  indicates  the  large  scale  vortical  structures 
produced  by  the  actuators  are  responsible  for  reattachment.  The  combination  of  momentum 
and  vorticity  injection  reattaches  the  how  over  the  airfoil  in  this  conhguration.  Spanwise 
vorticity  hux  observed  downstream  of  the  actuators  correlates  to  a  reattached  boundary 
layer.  It  should  be  observed  that  contour  plot  for  case  9C  shows  a  desirable  prohle  over 
the  entire  top  surface,  resulting  in  a  favorable  pressure  gradient  (attached  how)  achieved  by 
the  control  inputs  with  wall-normal  momentum  and  vorticity. 

Controlling  the  how  at  a  =  9°  requires  both  wall-normal  momentum  and  vorticity  injec- 
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C, 

Cv 


Surface  vorticity  flux  cry 
CTz 


Baseline 
0.5% 

0 


9C 
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Instantaneous  (Qj  Cp) 

Time- averaged  {ux  =  0,  Txy) 
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Figure  3.3:  For  each  case  at  a  =  9°,  the  instantaneous  flow  held  (top  left)  is  shown  with 
Q-criterion  colored  with  pressure.  The  time-averaged  zero-streamwise  velocity  {ux  =  0)  iso¬ 
surface  is  colored  with  spanwise  Reynolds  stress  Txy  (bottom  left).  Visualized  on  the  right 
are  the  wall- normal  vorticity  hux  components  in  the  streamwise  (top),  wall-normal  (middle), 
and  spanwise  (bottom)  directions. 
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tions  to  reattach  to  the  flow  based  on  the  representative  cases  considered  above.  Flow  control 
with  only  wall-normal  momentum  input  (case  9A)  perturbed  the  shear  layer,  but  does  not 
signiflcantly  impact  the  separated  wake.  The  addition  of  wall-normal  vorticity  reduces  the 
size  of  the  separated  flow  region.  From  case  9B  with  we.max/f^oo  =  0.95,  we  observe  a  decrease 
in  the  size  of  the  separated  flow  region.  Furthermore,  in  case  9C  with  M6),max/f^oo  =  1-26, 
separation  is  minimized  downstream  of  the  actuators  resulting  in  37%  drag  decrease  and 
31%  lift  increase.  With  all  the  cases  examined  herein,  we  summarize  the  influence  of  flow 
control  inputs  in  terms  of  the  aerodynamic  forces  experience  by  the  NACA  0012  airfoil  in 
the  following  sections. 
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Chapter  4 

Control  Input  Quantification 


Using  the  knowledge  gained  from  the  flow  fields  in  the  previous  section  let  us  further  elab¬ 
orate  on  the  flow  physics  and  quantify  the  momentum  introduced  to  the  boundary  layer 
as  a  function  of  actuator  input  parameters.  To  describe  the  control  input,  we  will  use  the 
wall-normal  and  azimuthal  velocities  input  by  the  actuators  as  well  as  the  geometry  of  the 
actuator.  We  separate  the  effects  of  momentum  directly  input  into  the  boundary  layer,  and 
momentum  introduced  to  the  boundary  layer  by  streamwise  vortices,  which  we  identify  as 
the  coefficient  of  momentum  and  coefficient  of  circulation.  Combining  these  two  coefficients 
results  in  what  we  will  refer  to  as  the  total  coefficient.  Using  the  following  quantihcation 
should  create  a  push  for  thorough  bench-top  characterization  of  different  actuators,  thus 
allowing  all  actuators  to  be  compared  in  a  similar  manner.  After  dehning  the  input  parame¬ 
ters,  we  examine  the  resulting  aerodynamic  forces  based  on  these  coefficients.  For  all  results 
presented,  force  induced  by  the  actuator  is  included  in  the  reported  forces. 


4.1  Coefficient  of  Momentum 

The  influence  of  momentum  directly  added  to  the  boundary  layer  is  quantihed  strictly  with 
the  wall-normal  velocity  component  and  actuator  radius.  By  convention,  we  use  the  mo¬ 
mentum  coefficient  which  is  the  ratio  between  the  momentum  input  by  the  actuator  to  the 
momentum  of  the  free  stream.  For  the  wall-normal  velocities  discussed  above,  the  momen¬ 
tum  added  was  not  adequate  to  reattach  the  separated  flow.  We  observed  in  the  following 
section,  increasing  the  wall-normal  velocity  eventually  reattaches  the  flow.  The  normal  ve¬ 
locity  used  in  this  study  result  in  coefficient  of  momentum  0(0.1%  —  1%),  which  is  of  similar 
magnitude  successfully  implemented  by  previous  studies  for  control  over  symmetric  airfoils 
(Deng  et  al,  2007;  Gilarranz  et  al,  2005;  Seifert  &  Pack,  1999;  You  et  al,  2008).  Next  the 
forces  are  plotted  with  versus  the  coefficient  of  momentum. 
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Figure  4.1:  The  coefficients  of  drag  (left)  and  lift  (right)  versus  the  coefficient  of  momentum. 
The  baseline  value  is  indicated  by - and  the  controlled  cases  are  pure  blowing,  Q;  co¬ 

rotating,  V,  and  counter-rotating,  A. 


4.1.1  Force  Coefficient 

The  resultant  drag  and  lift  forces  versus  the  coefficient  of  momentum  are  plotted  in  Fig.  4.1. 
The  information  contained  in  Fig.  4.1  is  as  follows;  hrst  we  distinguish  between  pure  blowing 
(O))  co-rotating  jets  (V),  and  counter-rotating  jets  (A).  Coefficient  of  momentum  is  only 
quantihed  by  wall-normal  velocity,  therefore  we  color  the  rotating  jet  cases  with  the  azimuthal 
velocity.  Additionally,  overlaid  on  these  plots  is  a  region  of  influence,  in  which  we  encompass 
the  resultant  forces  obtained  from  this  study.  The  region  of  influence  aids  in  understanding 
the  effect  of  the  coefficient  of  momentum  by  highlighting  the  range  of  forces  attained  for 
each  coefficient  of  momentum. 

When  flow  passes  over  the  airfoil  at  a  =  9°  and  Re  =  23,000,  the  drag  forces  show 
a  generally  decreasing  trend.  The  change  in  drag  force  is  minor  (AC'/j  ~  0.04),  but  the 
percentage  change  is  fairly  drastic  (AC/?  ~  35%).  While  almost  all  of  the  cases  reduced  drag 
on  the  airfoil,  lift  enhancement  is  achieved  with  appropriate  control  input.  For  coefficients 
of  momentum  0.15  <  C'^([%]  ^  0.45,  the  addition  of  wall-normal  vorticity  influences  the 
resulting  flow.  In  general,  for  a  constant  coefficient  of  momentum  (0.15%  <  C*^),  increasing 
ue  is  correlated  to  increasing  lift  and  decreasing  drag.  For  example,  cases  9A,  9B,  and  9C 
all  use  the  same  momentum  coefficient,  we  observe  that  increasing  azimuthal  velocity  or 
swirl  input  by  the  actuator  increases  lift  and  decreases  drag.  This  trend  is  noted  by  the 
diminishing  size  of  the  reverse  flow  region,  which  was  noted  when  comparing  cases  9 A,  9B, 
and  9C. 

From  these  hgures  we  ascertain  that  a  wall-normal  velocity  Un/U^o  >  0.625  effectively 
injects  the  wall-normal  vorticity  into  the  flow.  The  region  of  influence  illustrates  this  point, 
it  is  thin  for  small  momentum  coefficients  {un/Uoo  <  0.625).  It  appears  that  for  these  wall- 
normal  velocities,  the  swirl  added  to  the  boundary  layer  is  not  injected  near  the  edge  of 
boundary  layer  which  would  be  most  effective  to  mix  the  free-stream  and  boundary  layer. 
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The  gray  region  widens  with  increasing  coefficient  of  momentum,  and  now  the  vorticity  is 
being  injected  further  into  the  flow  by  the  wall-normal  momentum.  The  width  of  the  region 
of  influence  eventually  saturates  and  is  not  modified  by  the  larger  momentum  coefficients 
plotted  in  Fig.  4.1.  An  upper  limit  of  drag  decrease  and  lift  increase  appears  to  be  achieved 
for  cases  with  >  0.25,  and  wall-normal  vorticity.  Next,  let  us  define  the  coefficient  of 
circulation  and  evaluate  the  forces  as  a  function  of  the  coefficient  of  circulation. 


4.2  Coefficient  of  Circulation 


Considering  the  previous  discussion,  we  would  like  a  variable  similar  to  the  coefficient  of 
momentum  to  describe  the  circulation  input  to  the  flow.  In  the  previous  section  we  noted  that 
the  combination  of  wall-normal  momentum  and  swirl  reattached  the  flow.  The  increasing 
wall-normal  velocity  was  needed  to  advect  the  vorticity  further  into  the  boundary  layer. 
We  note  that  an  alternative  form  of  the  coefficient  of  momentum  is  the  dot  product  of  the 
non-dimensional  mass  flow  rate  (fnj/UooA)  and  the  blowing  ratio  (i?),  =  rUjUj/ 0.511^1^0. 

Therefore,  with  a  similar  thought  process,  we  multiply  the  non-dimensional  mass  flow  rate 
with  the  ratio  of  the  azimuthal  velocity  to  the  free-stream  velocity.  This  term  will  be  denoted 
as  the  coefficient  of  circulation  (Cr)  and  defined  as 


C'r  = 


rrijUe 


Un-n-rlue 

\UIA  ’ 


(4.1) 


and  similar  to  the  normal  velocity,  ug  =  ug^  max2/3.  While  this  might  not  appear  to  be  the 
obvious  choice  to  define  the  circulation  input  into  the  flow,  we  will  show  in  the  following 
that  it  works  well  to  describe  the  momentum  introduced  by  streamwise  vortices.  Coefficient 
of  circulation,  in  essence  is  similar  to  helicity  density,  which  is  the  dot  product  of  velocity 
and  vorticity.  Based  on  the  velocity  profiles,  the  circulation  input  by  the  actuator  is  T^  = 
QAiTUg^  max^"o/27.  Usiug  the  circulation  input  by  the  actuator  we  compute  the  non-dimensional 
vortex  strength  (8  <  T^/urh  <  53)  which  is  similar  to  vortex  strength  reported  by  Ashill 
et  al.  (2002)  (15  <  T/urh  <  45).  This  means  that  the  estimated  vortex  strength  values  in 
the  present  theoretical  study  are  realizable  with  present  vortex  generator  technology. 

In  the  previous  section,  we  identified  circulation  input  {f{ug))  and  penetration  depth 
{f{un/Uoo))  to  be  key  parameters  in  affecting  the  flow.  Here,  we  chose  to  non-dimensionalize 
the  circulation  input  such  that  r„  =  r„/f/oo4  =  327rM6iro/9f/oo4-  The  selection  to  non- 
dimensionalize  by  the  spanwise  spacing  of  the  actuators  is  because,  we  are  concerned  with 
streamwise  vortices  mixing  the  free-stream  and  boundary  layer.  A  streamwise  vortex  will 
impart  velocity  normal  to  the  boundary  layer  (mixing),  which  acts  between  the  two  actuators. 
Therefore  we  can  rewrite  the  coefficient  of  circulation  as 


C^r  = 


lU  c  ’ 
2 


(4.2) 


and  here  we  multiply  by  a  constant  C  to  account  for  the  fractional  term  in  circulation. 
The  remaining  terms  produce  the  blowing  ratio  {R  =  Un/Uoo)  which  has  been  identified 
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Figure  4.2:  The  coefficients  of  drag  (left)  and  lift  (right)  versus  the  coefficient  of  circulation. 
The  baseline  value  is  indicated  by - and  the  controlled  cases  are  pure  blowing,  Q,  co¬ 

rotating,  V,  and  counter-rotating,  A. 


previously,  and  when  multiplied  by  the  diameter  (ro/0.5)  we  obtain  a  coarse  estimate  of  the 
penetration  depth.  While  Muppidi  &  Mahesh  (2005)  establish  a  more  precise  estimation 
of  the  jet  height  {h)  based  on  the  downstream  distance  from  the  actuator  and  constants 
{A  and  5),  h/Rd  =  A{x/Rd)^,  we  note  the  reliance  on  the  blowing  ratio  and  diameter, 
Rd.  The  resulting  coefficient  of  circulation  is  the  non-dimensional  circulation  times  the  non- 
dimensional  height  of  the  actuator  (h  =  h/c).  Therefore  we  can  more  simply  express  the 
coefficient  of  circulation  as  Cr  =  CT h.  Next  we  will  assess  if  the  coefficient  of  circulation  is 
adequate  to  account  of  the  vorticity  input  by  the  actuator. 

4.2.1  Force  Coefficient 

Plotting  the  aerodynamic  forces  versus  4.2  (left)  and  we  note  some  general  trends;  lift  in¬ 
creases  and  drag  decreases  with  increasing  coefficient  of  circulation.  This  general  trend  does 
not  pertain  to  all  of  the  cases  but  does  largely  quantify  the  trend  observed.  For  example, 
cases  9C  and  9D  show  increased  lift  compared  to  lower  coefficients  of  circulation.  Although, 
further  increasing  the  coefficient  of  circulation  coincides  with  a  decrease  in  lift  force.  The 
aerodynamic  performance  parameter  {Cl/Cd)  is  plotted  in  Fig.  4.2  (right).  In  this  figure 
we  observe  that  there  is  a  drastic  difference  in  the  aerodynamic  performance,  it  ranges  from 
3.5  <  Cl/Cd  <  11.5.  From  this  plot  and  the  flow  field  results,  we  know  that  the  cases  with 
Cl/Cd  >  7  correlate  to  fully  attached  flow.  For  most  cases  it  is  apparent  based  on  the  forces 
the  result  of  the  flow  control  input.  With  the  above  quantification  of  swirl  input  to  the  flow, 
the  pure  blowing  cases  are  zero.  Therefore,  these  cases  are  on  the  y-a.xis  in  Fig.  4.2.  These 
plots  indicate  that  the  above  coefficient  does  begin  to  quantify  the  effects  of  vorticity  input, 
but  does  not  sufficiently  collapse  the  data  without  the  inclusion  of  wall-normal  momentum 
input.  Next,  we  include  the  influence  of  both  terms  to  account  for  all  of  the  momentum 
introduced  to  the  boundary  layer  by  the  control  input. 


0  9D  ^ 
9C  V 


11.0 

10.0 
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Figure  4.3:  The  drag  and  lift  forces  versus  the  total  coefficient.  The  baseline  value  is  indicated 
by - and  the  controlled  cases  are  pure  blowing,  Q;  co-rotating,  V,  and  counter-rotating, 

A. 


4.3  Coefficients  of  Momentum  and  Circulation 

The  coefficient  of  momentum  has  been  used  extensively  in  the  past  to  describe  the  influence 
of  actuators,  and  with  the  inclusion  of  the  coefficient  of  circulation  we  attempt  to  include 
another  component  that  affects  the  effectiveness  of  control  input.  Combining  these  two 
coefficients  (Ctotai  =  C^  +  Cr)  allows  us  to  more  accurately  quantify  the  total  momentum  in¬ 
troduced  to  the  boundary  layer,  which  can  aid  in  the  application  of  a  wide  range  of  actuators 
to  single  flow  control  scenario.  The  force  coefficients  versus  the  total  coefficient  are  plotted  in 
Fig.  4.3.  Additional  cases,  that  are  discussed  below,  have  been  added  to  this  plot  to  allow  for 
further  investigation  into  the  momentum  and  circulation  coefficient.  We  purposely  highlight 
the  pure  blowing  cases  with  hlled  circles,  because  these  values  are  identical  to  the  values  for 
pure  blowing  in  Fig.  4.1.  In  the  pure  blowing  cases  Cp  =  0,  and  Ctotai  =  C^,  therefore  we 
are  trying  to  use  the  coefficient  of  circulation  to  account  for  the  influence  of  vorticity  added 
to  flow;  shifting  the  values  to  align  more  closely  to  the  trend  produced  by  the  pure  blowing 
cases. 
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With  the  current  selection  of  +  Cy  we  obtain,  in  general,  two  prevailing  trends  exist 
which  correlate  to  the  separated  flow,  and  fully  reattached  flow.  Plotting  lift  versus  the 
total  coefficient,  three  regions  appear  in  Fig.  4.3.  For  smaller  control  input  (Ctotai  ^  0.6%), 
lift  force  decreases  for  the  majority  of  cases  which  is  due  adding  the  perturbation  to  the 
flow.  Though  the  perturbation  does  not  add  adequate  momentum  to  overcome  the  adverse 
pressure  gradient.  After  the  initial  drop  due  to  the  perturbation,  the  lift  forces  are  fairly 
consistent.  In  Fig.  4.4  this  region  is  very  flat  and  the  aerodynamic  performance  is  consistent 
for  the  pure  blowing  cases.  This  first  region  primarily  correlates  to  separated  flow,  and  is 
described  by  case  9A,  in  which  the  time-averaged  recirculation  region  is  approximately  the 
same  size  as  the  baseline  case.  We  also  note  that  in  this  region  there  are  cases  that  reattach 
the  flow. 

For  these  cases  which  are  encircled  in  Fig.  4.3,  flow  reattaches  when  the  majority  of  other 
cases  are  separated  with  this  control  input.  The  outlier  Case  9D,  is  plotted  in  Fig.  4.6  with 
the  same  visualizations  as  Fig.  3.3.  We  notice  that  this  case  eliminates  separation  similar 
to  case  9C.  The  zero-streamwise  velocity  iso-surface  shows  very  little  Reynolds  stress,  which 
indicates  that  the  flow  is  primarily  modified  by  the  large  structures  as  opposed  to  turbulent 
fluctuations.  The  signihcant  difference  for  the  circled  V  and  A  cases  is  the  blowing  ratio  is 
such  that  h{un/Uoo)/5  ^  1  based  on  the  relation  described  by  Muppidi  &  Mahesh  (2005). 
For  these  cases  that  reattach  the  flow  with  Ctotai  ^  0.6%,  the  interaction  of  wall-normal  and 
azimuthal  vorticity  are  ideal  to  induce  the  interaction  of  the  free  stream  and  boundary  layer, 
leading  to  reattached  flow.  Using  an  azimuthal  velocity  0=1  input  at  approximately  the 
boundary  layer  height  is  the  most  effective  location  for  the  streamwise  vortices  to  mix  the 
flow,  which  is  why  these  cases  reattach  the  flow  with  a  smaller  total  coefficient. 

Near  Ctotai  ~  0.65%  lift  begins  to  increase,  shown  in  Fig.  4.3,  and  correlates  to  the 
increase  in  aerodynamic  performance  in  Fig.  4.4.  In  the  vicinity  of  this  region  we  observe  a 
small  transition  region,  0.65  <  Ctotai  ^  0.8%,  which  corresponds  to  case  9B.  The  increase  in 
these  values  is  a  result  of  reducing  the  size  of  the  reverse  flow  region.  For  case  9B,  control 
input  is  impacting  the  separation  region,  but  flow  still  remains  separated  from  the  leading 
edge.  The  hnal  region  is  attached  flow,  Ctotai  ^  0.8%,  and  is  depicted  by  case  9C.  For  cases 
that  are  shown  in  this  region  we  can  begin  to  be  conhdent  that  any  device  selected  with 
Ctotai  ^  0.8%  will  reattach  the  flow.  All  of  the  cases  we  examined  with  a  total  coefficient 
greater  than  0.8  reattach  the  flow  and  result  in  similar  force  coefficients.  The  maximum 
attainable  increase  in  lift  and  drag  reduction  is  for  an  attached  flow,  so  further  increasing 
input  does  not  continue  increasing  the  aerodynamic  performance.  These  three  regions  are 
not  exact,  but  for  any  properly  quantified  actuators  would  indicate  how  much  input  is  needed 
to  reattach  the  flow. 

From  the  previous  cases  discussed  (9A,  9B,  9C,  9D)  and  our  initial  estimations  based  on 
C^,  we  did  not  attach  the  flow  with  pure  blowing.  Taking  the  knowledge  of  the  previous 
discussion  and  the  Ctotai  >  0.8  cut-off  into  consideration,  we  estimated  a  wall-normal  velocity 
that  would  reattach  the  flow  based  on  the  total  coefficient.  Using  this  value  (Ctotai  =  C^  = 
0.95%,  case  9E),  we  observe  in  Fig.  4.3  that  this  pure  blowing  case  follows  the  aforementioned 
trend.  As  expected  using  a  value  with  C^tai  >  0.8%  has  aerodynamic  forces  similar  to  cases 
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Figure  4.4:  The  aerodynamic  performance  versus  the  total  coefficient.  The  baseline  value 

is  indicated  by - and  the  controlled  cases  are  pure  blowing,  Qi  co-rotating,  V,  and 

counter-rotating,  A. 
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a 

Case 

'^n,max  / 

"^9  ^raax  /  U 

Rot.  dir. 

Cd 

Cl 

9° 

Baseline 

- 

- 

- 

0.118 

0.584 

9D 

0.945 

1.26 

CTR 

0.084 

0.673 

9E 

1.83 

0 

- 

0.084 

0.673 

9F  {IJc  =  0.2) 

1.78 

0 

- 

0.075 

0.725 

Table  4.1:  Select  control  parameters  for  flow  over  a  NACA0012  airfoil  at  a  =  9°.  The  case 
names  are  for  control  settings  that  are  used  in  the  text  herein. 


with  attached  flow.  The  resulting  flow  held  is  plotted  in  Fig.  4.6  (top)  and  we  observe  that 
the  how  is  indeed  attached  downstream  of  the  actuators.  The  increased  maximum  velocity 
increases  both  the  momentum  injected  into  the  how  leading  to  attached  how  for  this  case. 


4.4  Spanwise  Spacing 

The  total  coefficient  is  normalized  by  the  area  that  each  actuator  ahects,  thus  taking  into 
account  spanwise  spacing.  While  keeping  =  0.44%,  spacing  smaller  than  /^/c  <  0.1  was 
investigated.  The  resulting  input  values  are  0.15%  <  Ctotai  <  1-0%  were  used  in  an  attempt 
to  control  the  how  with  this  conhguration.  It  was  found  that  the  actuators  are  too  close 
together  and  interact  with  one  another.  This  interaction  causes  control  to  be  less  ehective 
and  no  cases  were  run  in  which  the  how  was  successfully  attached  for  l^jc  =  0.0667  and  0.05. 

Increasing  the  spacing  to  l^/c  =  0.2  allows  the  actuators  to  function  without  interfering 
with  the  adjacent  neighbor.  The  normal  velocity  is  selected  such  that  the  coefficient  of 
momentum  is  the  same  as  cases  run  for  l^/c  =  0.1.  To  keep  the  coefficient  of  momentum 
constant,  the  maximum  wall- normal  velocity  is  increased,  which  ohsets  the  increase  in  area. 
The  cases  run  at  l^jc  =  0.2  are  plotted  in  Fig.  4.5.  Drag  decreases  with  increasing  total 
coefficient,  while  lift  increases  once  how  is  attached.  With  the  cases  run,  we  observe  that  the 
data  adheres  to  the  previously  reported  trend.  For  this  spacing  how  reattaches  at  a  lower 
total  coefficient.  We  observe  in  Fig.  4.6  how  is  attached  downstream  of  the  actuator  for  case 
9F,  which  has  how  features  similar  to  9E.  This  indicates  that  the  spacing  of  the  actuators 
an  important  parameter  when  attempting  to  optimize  how  control  input. 
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Table  4.2:  Time  history  of  the  drag  (top  left)  and  lift  (bottom  left)  coefficient  for  a  pure 

blowing  actuator  at  ^jc  =  0.05  The  baseline  force  is  shown  with  .  .  The  flow  held 

images  correspond  to  the  •  in  the  time  history  plots,  and  from  top  to  bottom  the  images 
are  sequential  in  time. 
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Figure  4.5:  The  drag  and  lift  forces  (left)  and  aerodynamic  performance  (right)  versus  the 

total  coefficient.  The  baseline  value  is  indicated  by - and  the  controlled  cases  are  pure 

blowing,  O)  co-rotating,  V,  and  counter-rotating,  A. 
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Figure  4.6:  For  each  control  case  in  Table  4.1  the  instantaneous  flow  fleld  (top  left)  is 
shown  with  Q-criterion  colored  with  pressure.  The  time-averaged  zero-streamwise  velocity 
iso-surface  is  colored  with  spanwise  Reynolds  stress  (bottom  left).  The  wall  normal  vorticity 
flux  is  to  the  right  from  top  to  bottom  is  the  streamwise,  wall-normal,  and  spanwise  vorticity. 
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Chapter  5 

Global  Stability  Analysis 


Understanding  the  temporal  and  spatial  evolution  of  perturbations  can  aid  in  the  placement 
of  flow  control  devices  as  well  as  further  the  understanding  the  effects  perturbations  have 
on  the  overall  flow  held.  Linear  theory  predicates  that  how  helds  can  be  comprised  of 
stable  (exponential  decay),  unstable  (exponential  growth),  and/or  neutrally  stable  (remains 
constant)  modes;  identihcation  of  these  modes  can  be  done  with  several  diherent  approaches. 
These  approaches  include  dynamic  mode  decomposition  (Rowley  et  al,  2009;  Schmid,  2010), 
initial  value  problems  (Brehm  &  Fasel,  2011),  and  global  stability  analysis  (Theohlis,  2011, 
2003),  each  with  advantages  and  limitations.  For  example,  while  the  initial  value  problems 
can  be  performed  at  a  lower  computational  cost,  yet  the  ability  to  resolve  the  eigenspectrum 
is  greatly  reduced.  Global  stability  analysis  provides  a  broader  range  of  eigenvalues  at  an 
increased  computational  cost. 

5.1  Bi-Global  Stability  Formulation 

Performing  linear  stability  analysis  hrst  begins  with  the  incompressible  Navier-Stokes  equa¬ 
tions. 


dU  „  „  1  , 

—  +  u  ■  Vu  = -Vp  + —V^u,  (5.1) 

ot  Re 

V  ■u  =  0,  (5.2) 

in  which  we  use  the  non-dimensional  form,  noted  by  the  Reynolds  number.  Re  =  Uool/voo- 
Equations  5.1  and  5.2  are  linearized  such  that  the  solution  {u  and  p)  is  the  sum  of  a  base 
state  {u  and  p)  and  perturbation  {u'  and  p')  (Schmid,  2007;  Theohlis,  2011,  2003).  In  the 
linearization  process,  two  assumptions  are  made:  the  base  state  {u,  p)  satishes  the  Navier- 
Stokes  equations,  and  the  perturbations  {u',p')  are  relatively  small,  i.e.,  |it'|/|lZ|  =  0{e)  1. 

When  u  +  u'  and  p  +  p'  are  substituted  into  Eqs.  5.1  and  5.2,  and  higher  order  terms  (G(e^)) 


31 
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are  removed,  we  arrive  at  the  linearized  Navier-Stokes  equations 


dt  Re 


(5.3) 

(5.4) 


V  •  ti'  =  0. 


From  the  above  equations,  we  not  that  the  advection  term  provides  the  only  mechanism  for 


the  base  flow  to  interact  with  the  perturbation.  Once  perturbations  grow  large  (|w'|/|w|  >  e). 


the  linearized  equations  are  no  longer  valid  due  to  the  sizeable  nonlinear  interactions  amongst 
the  perturbations. 

In  order  to  simplify  the  problem,  assumptions  are  made  about  spatial  characteristics  of 
the  flow.  The  simplihcation  assumes  a  wave-like  propagation  in  one  or  two  spatial  directions. 
For  flows  varying  in  two  spatial  directions,  assuming  a  wave-like  pattern  in  the  third  is 
considered  bi-global  stability.  To  perform  bi-global  stability  analysis,  the  perturbation  q' 
=  [u',p']'^  is  assumed  to  have  a  solution  of  the  form 


q'(x,  y,  z,  t)  =  q{x,  y)  exp[!(/32  -  uit)]  +  c.c., 


(5.5) 


and  the  perturbation  is  assumed  to  be  wavelike  in  the  spanwise  direction  with  a  spanwise 
wavenumber  /3  =  27r jL^,  and  is  the  spanwise  wavelength.  Temporally,  u  =  Ur  +  ioJi 
is  the  complex  eigenvalue  for  the  two-dimensional  eigenmode  q.  Based  on  Eq.  5.5,  the 
eigenvalues  contain  information  about  the  exponential  growth  and  decay  rates  (cuj),  as  well 
as  the  frequency  of  the  global  modes  {oJr)-  The  growth  {ui  >  0)  or  decay  {ui  <  0)  rate 
is  of  particular  interest  because,  as  mentioned  previously,  flow  control  can  effectively  take 
advantage  of  inherent  instabilities  in  the  flow.  Herein,  we  refer  to  modes  with  a  positive 
growth  rate  as  unstable,  and  negative  as  stable  (wj  =  0,  is  neutrally  stable).  Additionally, 
modifying  the  frequency  of  perturbations  can  directly  influence  the  growth  and  decay  rate 
of  specific  modes.  Based  on  the  frequency  content  of  modes,  we  classify  a  mode  as  steady  if 
cUr  =  0  or  unsteady  if  Ur  7^  0.  Unsteady  modes  come  in  complex  conjugates  pairs. 

5.1.1  Eigenvalue  Problem 

To  perform  an  eigenanalysis  of  a  fluid  flow,  we  seek  to  derive  the  algebraic  form  of  the  lin¬ 
earized  Navier-Stokes  equations.  Substituting  the  solution  form  (Eq.  5.5)  into  the  linearized 
Navier-Stokes  equations  (Eqs.  5.3  and  5.4)  results  in 


—lUV, 


—lUU, 


—lUW, 


(5.6) 


(5.7) 


(5.8) 


du  dv 
dx  dy 


+  ij3w  =  0. 


(5.9) 
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We  can  formulate  these  equations  as  Aq  =  —iXBq,  where  A  is  the  linear  operator  (left-hand 
side  of  Eqs.  5.6  -  5.9),  and  A  is  the  eigenvalue  to  be  solved  (a;).  Matrix  B  is  of  the  form 


B 


A  0  ■■■O' 

0  1  ; 

;  1  0 

0  ■■■  0  0 


(5.10) 


which  by  simple  inspection,  is  not  invertible  {det{B)  =  0).  Thus,  the  shift  and  invert  strategy 
is  utilized,  and  the  eigenvalue  problem  becomes  {A  —  aB)~^Bq  =  Xq  (Natarajan  &  Acrivos, 
1993).  The  new  problem  solves  different  eigenvalues  (A),  which  are  related  to  the  original 
problem  by  A  =  1/(A  —  a).  This  particular  algorithm  searches  for  eigenvalues  closest  to  the 
shift  parameter  (a).  The  shift  value  can  be  modified  to  obtain  results  for  different  regions 
of  the  eigenspectrum. 

We  take  advantage  of  the  domain  decomposition  and  discrete  operators  from  the  Charles 
software  package.  As  mentioned  previously,  the  finite  volume  operators  are  second  order 
accurate  (Ham  &  laccarino,  2004;  Ham  et  ai,  2006).  For  this  discretization  scheme,  the  non¬ 
zero  values  of  the  matrix  are  illustrated  in  Fig.  5.1.  In  the  present  study,  a  two-dimensional 
base  state  will  result  in  a  non-zero  value  in  only  two  directions  (u,  v  ^  0,  w  =  0  ),  and 
the  matrix  structure  is  plotted  Fig.  5.1.  Due  to  the  size  and  sparseness  of  the  matrices 
A  G  and  B  G  ,  a  matrix-free  method  to  solve  the  eigenvalue  problem  is  desired 

to  reduce  memory  requirements  for  the  large-scale  problem.  The  two-dimensional  mesh  for 
the  NACA  0012  airfoil  has  (9(10^)  points,  and  the  size  of  the  problem  N  =  dxnumber  of 
grid  points.  The  reverse  communication  interface  of  ARPACK  (Lehoucq  et  ai,  1998)  is  used 
to  accommodate  the  matrix-free  Arnoldi  algorithm.  An  iterative  solver  is  used  to  invert  the 
large  sparse  matrix  {A  —  (tB)~^\  Bi-Conjugate  Gradient  Stabilized  method  (Bi-CGSTAB) 
(Barrett  et  ai,  1993)  is  implemented  for  the  matrix  inversion. 

The  implementation  of  the  aforementioned  procedure  can  be  simplihed  and  explained  in 
Algorithm  1.  The  procedure  hinges  primarily  on  the  ARPACK  subroutine  dnaupd,  which  is 
the  serial  implicitly  restarted  Arnoldi  algorithm  (Lehoucq  et  ai,  1998).  The  parallel  version 
is  pdnaupd.  The  implicitly  restarted  Arnoldi  algorithm,  creates  an  output  {b)  and  requests 
an  input  (x)  based  on  the  operator  {A).  As  mentioned  previously,  the  returned  values  (a;), 
are  solved  with  the  Bi-CGSTAB.  The  convergence  to  the  proper  eigenvalues  is  dependent  on 
the  convergence  criteria  of  the  two  iterative  processes.  According  to  the  documentation  of 
the  ARPACK  subroutine  used  (Lehoucq  et  ai,  1998),  the  Bi-CGSTAB  convergence  tolerance 
should  be  less  than  the  tolerance  of  the  Arnoldi  method  (to/Bi-cc  <  to/Amoidi)- 
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Figure  5.1:  Structure  of  the  full  matrix  form  of  Eqs.  5.6  -  5.9  on  a  16  x  16  domain  with  peri¬ 
odic  boundary  conditions.  The  colors  indicate  the  non-zero  real  and  imaginary  components 
of  real  and  imaginary  components  are  represented  by  ■  ,  and  ■  ,  respectively. 


Algorithm  1  Matrix  Free  Bi-Global  Stability  Analysis 

1:  while  errAmoidi  >  tol Amoidi  do 

2:  Call  Arnold!  algorithm  (znaupd/pznaupd) 

3:  Interpret  output  from  Arnold!  algorithm 

4:  Set-up  the  problem  in  the  form  Ax  =  b 

5:  Provide  an  initial  guess  for  xq 

6:  Compute  initial  residual 

7:  while  errBi-cGSTAB  >  tolBi-CG  do 

8:  Update  initial  guess  based  on  residual  (xj_i  x*) 

9:  Compute  new  residual  based  on  new  solution 

10:  Update  error  (errsi-cG  =  \\Axi  —  5||2) 

11:  end  while 

12:  end  while 

13:  Call  Arnold!  post-processing  routing  (zneupd/pzneupd) 


5.1.2  Boundary  Conditions 

Traditional  boundary  conditions  are  enforced,  but  special  consideration  for  stability  analysis 
has  to  be  taken  into  account  for  both  the  base  flow  and  perturbation.  Dirichlet  or  Neumann 
boundary  conditions  can  be  implemented  for  velocity  and  pressure  perturbations. 
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Wall 


For  no-slip  boundaries,  velocity  at  the  wall  is  zero,  U bc  =  ^  =  u  =  The  base  flow 
should  satisfy  the  boundary  condition,  so  the  perturbation  velocity  at  the  wall  is  Dirichlet 
and  specified  to  be  h  =  h  =  ti)  =  0.  The  compatibility  pressure  boundary  conditions  are 
specified  to  be  Neumann, 

|  =  (5,11) 

which  is  reported  to  be  successful  in  the  previous  study  by  Theofllis  et  al.  (2004). 

Inflow 

Depending  on  the  structure  of  the  domain,  an  inflow  condition  may  be  more  appropriate  or 
applicable  as  opposed  to  the  far  held  condition.  In  the  present  cases,  the  inflow  is  laminar 
{Uoo  =  1),  and  with  that  knowledge  of  the  flow  field,  perturbations  should  approach  zero 
at  the  inlet  of  the  domain.  Therefore,  this  additional  boundary  condition  was  created  for 
perturbation  velocity  and  pressure  to  be  zero  {q'  =  0). 


Outflow  and  Far  Field 


Downstream  of  the  airfoil,  an  appropriate  boundary  condition  is  necessary  to  allow  pertur¬ 
bations  to  leave  the  domain.  Similarly,  approaching  the  far  held,  exact  boundary  conditions 
are  unknown,  therefore  the  outflow  boundary  condition  is  compatible  for  boundaries  consid¬ 
ered  as  far  field.  In  open  flow,  perturbations  are  advected  downstream  and  will  impact  the 
boundary  of  the  finite  computational  domain.  Boundary  information  is  obtained  from  the 
interior,  by  means  of  linear  extrapolation.  Values  of  the  interior  perturbation  and  derivative 
are  used  to  calculate  the  boundary  values.  To  accommodate  unstructured  domains,  multiple 
interior  nodes  (n  is  the  number  of  nodes)  are  used  to  extrapolate  the  solution  for  a  single 
boundary  node,  which  is  done  as  follows 


1 

ri  1 

i=l  Axi 


(5.12) 


Weighting  the  interior  nodes  to  the  inverse  of  distance  allows  for  nodes  closer  to  the  boundary 
node  to  be  more  influential.  The  distance  between  the  inner  node  and  boundary  point  are 
defined  as  Ax^  =  x^,  —  Xi.  Linear  extrapolation  at  the  boundary  has  worked  in  the  past  for 
multiple  studies  and  was  found  to  be  a  robust  approximation  for  boundaries  in  open  flow 
(Kitsios  et  al,  2009;  Rodriguez,  2010). 


5.1.3  Base  flow 

In  the  formulation  of  stability  analysis,  a  base  state  is  required  for  the  linear  Navier-Stokes 
equations.  To  solve  for  the  base  state,  direct  numerical  simulations  are  performed  using 
the  finite  volume  solver  Charles  (Cliff  for  Incompressible  flow).  Boundary  conditions  used 
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are  slip,  no-slip,  inlet,  and  outlet.  For  steady  flow,  solving  for  the  base  flow  is  trivial.  A 
simulation  can  be  performed  allowing  transients  to  decay,  resulting  in  the  final  state  being 
equivalent  to  the  stable  base  state.  It  is  more  difficult  to  obtain  the  base  state  for  an  unstable 
flow,  and  we  outline  the  formulation  used  below. 

Unstable  steady  state 

In  the  present  study,  the  primary  intent  is  to  perform  stability  analysis  with  respect  to  the 
unstable  steady  state.  Note  that  the  base  flow  is  in  general  not  equal  to  the  time-averaged 
flow  field.  Selective  frequency  damping  is  used  to  find  the  unstable  steady  state  and  can  be 
represented  as 

rjqt 

—  =  NS{u)-x{u-u),  (5.13) 

where  N S  is  the  right-hand  side  of  the  Navier-Stokes  operator,  y  is  the  gain,  and  if  is  a  time- 
filtered,  nonlinear  flow  field  (Akervik  et  al,  2006).  The  forcing  term  added  to  the  right  hand 
side  of  the  Navier-Stokes  equations  forces  the  flow  to  the  unstable  steady  state.  A  second- 
order  low-pass  Butterworth  filter  is  implemented  to  allow  the  filtered  solution  to  approach 
the  base  state.  The  cut-off  frequency  is  selected  to  be  below  the  dominant  (characteristic) 
frequency  present  in  the  flow.  This  case  is  run  until  max{du/dt)  <  tol,  where  tol  is  a 
tolerance  selected  by  the  user. 

5.2  BiGlobal  Stability  Analysis 

We  validate  the  above  formulation  using  the  following  two  problems:  lid-driven  cavity  and 
separated  flow  over  a  NACA  0015  at  a  =  18°.  First,  we  quantitatively  verify  the  results  by 
comparing  the  least  stable  eigenvalues  to  previous  lid-driven  cavity  studies.  Qualitatively, 
we  compare  the  eigenmodes  of  separated  flow  over  an  airfoil  to  verify  the  results  for  open 
flow  cases. 

5.2.1  Lid-Driven  Cavity 

BiGlobal  stability  analysis  of  lid-driven  cavity  flow  has  been  investigated  thoroughly  in  past 
studies.  The  steady  solution  {u{t  — )■  oo))  of  the  non-linear  direct  numerical  simulation  of  lid- 
driven  cavity  is  used  for  the  base  flow.  The  base  flow  for  Re  =  100  is  plotted  in  Fig.  5.2  (left). 
For  validation,  the  leading  eigenvalue (s)  for  two  different  spanwise  wavelengths  are  compared 
with  results  reported  by  Theofilis  et  al.  (2004),  see  Table  5.1.  Changing  the  spanwise  wave 
length  (/3),  modifies  the  characteristics  of  the  leading  eigenvalue.  The  leading  eigenvalue  for 
the  both  spanwise  wavelengths  is  stable  (aj  <  0),  and  is  steady  (cu^  =  0)  for  /3  =  1  and 
unsteady  {ui  ^  0)  for  (3  =  5.  The  first  96  eigenvalues  for  Re  =  200  and  (3  =  1  are  shown  in 
Fig.  5.2  (right). 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Table  5.1:  Most  unstable  eigenvalue  for  lid-driven  cavity  flow  at  Re  =  200. 


(5  =  1 

LO 

II 

CQ, 

Theofilis  et  al.  (2004) 

0.000  -  0.3297i 

±0.4260  -  0.3404i 

Present  Results 

0.000  -  0.3298i 

±0.4274  -  0.3410i 

Figure  5.2:  Base  flow  and  96  closest  to  the  shift  value  a  =  1.  The  Reynolds  number  and 
spanwise  wave  number  are  Re  =  200  and  (3  =  1,  respectively. 
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5.2.2  NACA  0015,  a  =  18° 


A  NACA  0015  airfoil  at  a  =  18°  is  investigated  at  low  Reynolds  numbers  (i?eC>(100))  by 
Kitsios  et  al.  (2009)  and  we  will  use  this  set-up  to  verify  the  correct  implementation  of  the 
code  for  open  flow  scenarios.  For  comparison  at  higher  Reynolds  numbers  {ReO{10^))  we 
compare  the  results  of  the  present  solver  to  those  of  Zhang  &  Samtaney  (2016).  Though 
the  study  by  Zhang  &  Samtaney  (2016)  is  with  a  different  airfoil  (NACA  0012)  and  angle 
of  attack  («  =  16°),  the  base  state  for  the  separated  flows  is  similar.  Thus,  we  qualitatively 
compare  the  results  to  ensure  the  solver  is  capable  of  handling  open  flow  and  increased 
Reynolds  number.  We  compare  the  results  of  flow  at  Re  =  200,  600,  and  1,  000  and  (3  =  1. 

Base  Flow 

For  the  following  cases,  the  base  flow  is  the  unstable  steady  state,  except  for  Re  =  200 
which  we  will  notice  is  inherently  steady.  The  gain  used  in  all  of  the  cases  using  selective 
frequency  damping  is  y  =  1.0.  The  size  of  the  two-dimensional  computational  domain 
{x/c,y/c)  G  [—20,25]  x  [—20,20].  The  boundary  conditions  at  the  inlet  is  uniform,  laminar 
flow,  Uco  =  [1,  0].  Top  and  bottom  are  slip  boundaries.  The  airfoil  near  the  center  of  the 
domain  is  set  to  be  a  wall,  i.e.  no-slip  and  no-penetration  boundary.  The  outflow  boundary 
is  a  convective  outlet,  which  allows  structures  to  leave  the  domain  without  affecting  the  near 
held  solution.  The  outhow  boundary  is  described  by  ^  17c||  =  0  and  the  convective 

velocity  (17c)  is  dehned  as  the  mean  velocity  at  the  outlet. 

Computing  the  base  how  at  Re  =  200,  we  observe  a  stable,  separated  how.  In  Fig.  5.3 
(top),  how  separates  from  the  leading  edge  and  produces  a  primary  recirculation  region. 
From  the  trailing  edge  a  smaller  secondary  recirculation  region  is  generated.  The  center  of 
the  primary  recirculation  region  is  relatively  close  to  the  airfoil.  At  Re  =  600,  the  increase 
in  Reynolds  increases  the  size  of  both  the  recirculation  regions  behind  the  airfoil.  From 
Fig.  5.3  we  observe  further  increasing  the  Reynolds  number  gradually  increases  the  size  of 
the  recirculation  regions.  The  distance  of  the  two  recirculation  regions  are  a  similar  distance 
from  the  airfoil. 

Eigenmodes 

Using  the  base  hows  described  above,  we  compute  the  eigenvalues  and  eigenvectors  for  the 
BiGlobal  stability  formulation.  At  Re  =  200,  the  how  is  steady,  and  as  a  result  there  are 
no  unstable  eigenvalues.  In  Fig.  5.3,  the  most  unstable  modes  are  shown  and  at  Re  =  200; 
we  note  that  the  mode  has  a  negative  growth  rate.  The  colored  contour  visualizes  the  x- 
direction  perturbation  velocity,  which  is  similar  to  the  wake  mode  observed  in  separated  how, 
such  as  circular  cylinder  (Noack  &  Eckelmann,  1994).  Increasing  Reynolds  number,  both 
Re  =  600  and  1,  000  have  unstable  modes,  whose  eigenvectors  are  similar  that  of  Re  =  200. 
The  increased  Reynolds  number  produces  cases  that  begin  to  shed,  which  correlates  to 
the  emergence  of  unstable  modes.  At  Re  =  600,  there  are  two  unstable  modes:  a  steady 
mode  and  unsteady  mode  (plotted  in  Fig.  5.3).  Further  increasing  the  Reynolds  number 
{Re  =  1000),  an  additional  wake  mode  is  excited,  resulting  in  a  near  and  far  wake  mode. 
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Re  =  200,  u  =  ±2.34  -  iO.295 


Re  =  1000,  a;  =  ±2.23  ±  *0.271 


Figure  5.3:  Streamlines  of  the  base  flow  for  each  case  Reynolds  number,  accompanied  by  the 
contour  of  Ux- 


The  two  wake  modes  and  steady  mode  are  consistent  with  the  observations  of  Zhang  & 
Samtaney  (2016)  at  the  same  Reynolds  number.  The  agreement  between  the  two  cases 
shows  that  the  formulation  of  matrix-free,  bi-global  stability  eigenvalue  problem  appears 
to  be  correct.  With  direct  numerical  simulation  these  were  the  largest  Reynolds  number 
cases  that  were  attained  in  this  study  using  biglobal  stability  analysis.  Further  increasing 
Reynolds  number  requires  drastically  increased  computational  resources. 
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Chapter  6 

Network  Analysis 


6.1  Network-theoretic  framework 


To  assess  and  describe  the  vortical  interactions  in  the  flow  fleld,  we  utilize  a  weighted  network 
(graph).  The  deflnition  of  a  network  (graph)  Q  requires  sets  of  vertices  (nodes)  V,  edges  £, 
and  weights  W  (Newman,  2010).  With  these  three  components  defined,  a  graph  can  be 
uniquely  determined,  i.e.,  Q  =  ^(V,£^,  W).  The  nodes  V  in  this  study  are  taken  to  be  the 
vortical  elements  and  the  edges  £  represent  the  vortical  interactions  between  those  vortical 
elements.  The  edge  weights  W  quantify  the  strengths  of  the  vortical  interactions.  Given  n 
nodes,  a  collection  of  the  weights  Wij  in  the  form  of  a  matrix  A  G  with 


A,, 


u!ij  A{i,j)e£ 
0  otherwise. 


(6.1) 


is  called  the  adjacency  matrix  and  is  used  to  describe  the  network  connectivity.  In  the  above 
deflnition,  Aij  is  set  to  the  edge  weight  Wij  if  there  exists  an  edge  (interaction)  between  nodes 
i  and  j.  Details  on  the  fundamental  concepts  involved  in  network  theory  can  be  found  in 
Newman  (2010)  and  Dorogovtsev  (2010)  with  descriptions  of  vortical-interaction  networks  in 
Nair  &  Taira  (2015).  The  diagonal  entries  of  the  adjacency  matrix  relate  to  the  weight  of  the 
self-connecting  (loops)  edges.  The  adjacency  matrix  of  an  undirected  graph  is  symmetric. 
The  degree  ki  of  a  vertex  i  represents  the  summation  of  the  weights  of  the  edges  connected 
to  it  given  by  =  '^f=i[Ag]ij.  Another  important  matrix  in  graph  theory  is  the  graph 
Laplacian  matrix  Lg  G  which  is  given  by 


[Lg]ij 


{ki  a  {i,j)  e  E  and  {i  =  j) 

-Wij  if  {i,j)  G  E  and  {i  7^  j) 
0  otherwise. 


(6.2) 


The  graph  Laplacian  matrix  can  also  be  deduced  from  the  adjacency  matrix,  Lg  =  Dg  —  Ag, 
where  Dg  G  is  a  diagonal  matrix  with  elements  equal  to  degrees  of  vertices,  Dg  = 

diag([fci](^^).  For  undirected  networks,  it  is  symmetric  and  positive  semi  definite  (singular). 


40 
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Figure  6.1;  (left)  An  example  weighted  graph;  and  (right)  a  complete  graph  /C5. 


This  Laplacian  matrix  is  a  discrete  analog  of  the  negative  continuous  Laplacian  operator 
(— V^)  and  is  also  called  as  the  discrete  Laplacian.  It  is  naturally  dehned  by  its  quadratic 
form.  If  we  have  a  linear  system  of  the  form  Lgx  =  b,  for  the  vector  x  G  the  Laplacian 
quadratic  form  of  a  weighted  graph  Q  is  given  by  Mohar  (1991) 


x'^Lgx  = 


(6,3) 


The  discrete  Laplacian  is  a  smoothness  indicator  of  x  over  the  edges  in  Q.  The  Laplacian 
quadratic  form  becomes  large  as  x  jumps  over  the  edges  of  Q.  The  dehnitions  of  adjacency 
matrix  and  graph  Laplacian  form  the  building  blocks  of  the  graph-theoretic  framework.  The 
Laplacian  quadratic  form  will  be  used  below  to  introduce  the  notion  of  spectral  similarity  of 
graphs  and  graph  sparsihcation. 

An  example  of  a  weighted  graph  is  shown  in  Fig.  6.1  (left)  with  weights  associated  with 
the  edges  displayed.  A  complete  graph  /C^v  with  N  =  5,  shown  in  Fig.  6.1  (right),  has  all  the 
vertices  connected  to  each  other.  In  other  words,  this  graph  has  N  vertices  with  complete  set 
of  possible  N{N  —  l)/2  edges  without  any  self-loops.  The  adjacency  matrix  for  an  example 
graph  Ag 


and  the  complete  graph  Aycg 

shown  in  Fig.  6.1  are  given  by 
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1 
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/ 

where  ^  is  a  weighted  graph  and  /C5  is  an  unweighted  complete  graph. 
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Vortical  interactions  on  a  graph 

To  extract  the  network  structure  of  the  flow,  we  quantify  the  interactions  between  fluid 
elements  based  on  the  vortical  interactions.  The  velocity  u  at  position  x  induced  by  the 
vorticity  distribution  u:  of  the  flow  is 


u{x,  t) 


1 

47r 


u>(x,  t)  X  (x  —  x) 

- ^ - VT5 - dx. 

\x  —  a;  h 


(6,4) 


A  vortex  is  advected  by  the  induced  velocity  from  the  other  vortices  and  not  by  its  own 
velocity.  Each  vortex  is  influenced  by  all  other  vortices,  which  means  that  each  and  every 
vortex  has  a  connection  to  all  other  vortices  without  any  self-loops.  Thus,  the  discrete  system 
of  N  point  vortices  can  be  represented  by  a  weighted  complete  graph,  /Cat.  The  vertices  or 
nodes  of  the  graph  represent  the  discrete  point  vortices  and  the  edge  weights  represent 
the  strength  of  the  connections  between  them.  The  motion  of  point  vortices  is  influenced 
by  the  strengths  (circulation)  of  the  individual  point  vortices  and  the  distances  between 
them.  The  shorter  the  distance  between  the  point  vortices,  the  higher  is  the  influence  of 
the  vortices  on  each  other.  Also,  a  vortex  with  higher  strength  has  more  ability  to  influence 
the  surrounding  vortices  compared  to  that  with  lower  strength.  We  assign  the  weight  for 
the  interaction  between  two  vortices  to  be  dependent  on  the  characteristic  vortex  induced 
velocity  based  on  the  strengths  of  the  two  vortices  and  the  distance  between  them.  Thus, 
this  implies  that  the  edge  weight  for  the  vortex  network  should  be  proportional  to  n/r,  where 
K  is  the  characteristic  strength  of  the  two  vortices  and  r  is  the  distance  between  them  within 
the  network.  The  magnitude  of  the  induced  velocity  from  fluid  element  i  on  another  element 
j  reduces  from  Eq.  (6.4)  to 


Ui^j 


\  I 

2'K\Xi  —  Xj 


(6.5) 


Based  on  Eq.  (6.5),  we  dehne  the  network  adjacency  matrix  as  the  average  induced 
velocity 

h{ui^^  +  Uj^i)  \ii  ^  j 
0  otherwise 


A-.  =  ) 


(6.6) 


to  quantify  the  magnitude  of  interaction  between  fluid  elements  i  and  j  (Nair  &  Taira, 
2015).  Alternatively,  the  geometric  mean  of  induced  velocity  can  be  used  to  dehne  the 
network  adjacency  matrix. 


[Ag]  ij 


if  {i,j)  G  E  and  i  ^  j 
otherwise. 


(6.7) 


Note  that  an  element  cannot  impose  velocity  upon  itself,  which  is  captured  by  the  null  entry 
along  the  diagonal  of  the  adjacency  matrix. 
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6.2  Graph  sparsification 

In  order  to  identify  key  vortical  interactions  involved  in  nnsteady  flnid  flows,  we  construct 
sparse  approximation  of  dense  vortical  graphs.  Graph  sparsihcation  is  useful  for  design¬ 
ing  computationally  efhcient  algorithms  and  helps  in  identifying  representative  edges  and 
associated  weights  (Spielman  &  Teng,  2011).  Graph  similarity,  which  forms  the  basis  for 
graph  sparsihcation,  can  be  derived  in  a  number  of  ways  depending  on  the  desired  similar¬ 
ity  properties.  Distance  similarity  of  graphs  (Peleg  &  Ullnian,  1989)  can  be  achieved  by 
sparse  graphs  called  spanners  that  keep  the  same  shortest-path  distance  between  each  pair 
of  vertices  as  in  the  original  graph  while  cut  similarity  (Bencziir  &  Karger,  1996)  can  be 
achieved  by  maintaining  the  weight  of  the  edges  leaving  a  set  of  vertices  of  the  sparse  graph 
to  be  the  same  as  that  of  the  original  graph.  A  much  stronger  notion  of  similarity  referred 
to  as  spectral  similarity  was  introduced  by  Spielman  &  Teng  (2011).  Spectral  similarity  is 
closely  tied  to  the  Laplacian  quadratic  form  of  the  graphs  as  dehned  by  equation  (6.3).  The 
concept  of  spectral  similarity  directly  leads  to  spectral  sparsihcation  of  graphs;  that  is  to 
create  sparse  graphs,  which  are  spectrally  similar  to  the  original  graph. 

Spectral  sparsihcation  is  a  more  general  abstraction  than  cut  sparsihers  and  maintains 
spectral  similarity  between  the  sparsihed  and  original  graphs.  In  particular,  spectral  spar¬ 
sihcation  can  remove  some  of  the  edges  in  the  graph,  while  maintaining  similar  adjacency 
and  Laplacian  eigenspectra.  One  of  the  key  features  of  spectral  sparsihcation  is  that  it  keeps 
the  sum  of  the  weights  leaving  the  vertex  of  a  graph  constant.  A  spectral  sparsiher  is  a 
subgraph  of  the  original  graph  whose  Laplacian  quadratic  form  is  approximately  the  same 
as  the  original  graph  (Spielman  &  Teng,  2011). 

Sparsihcation  involves  the  creation  of  a  sparse  graph  Qs  from  the  original  graph  Q  based 
on  an  approximation  order  of  e.  The  quadratic  form  induced  by  graph  Laplacian  of  Q  is 
maintained  upto  a  multiplicative  (1  ±  e)  factor  by  spectral  sparsihcation  (Kelner  &  Levin, 
2011).  Thus,  the  sparse  graph  Qs  is  a  (l±e)-spectral  approximation  of  Q.  The  approximation 
order  e  can  vary  from  zero  to  unity.  The  approximation  with  e  =  0  indicates  that  Qs  is  same 
as  original  graph  Q  and  none  of  the  edges  are  sparsihed,  while  e  =  1  relaxes  the  quadratic 
form  induced  by  the  sparse  graph  to  within  twice  of  that  induced  by  the  original  graph.  An 
approximation  of  e  =  1  leads  to  a  heavily  sparsihed  graph.  Denoting  the  Laplacian  matrices 
of  Q  and  Qs  by  Lg  and  respectively,  the  spectrally  sparsihed  Laplacian  satishes 

(1  —  e)x^Lgx  <  x^LggX  <  (1  +  e)x'^Lgx  (6.8) 

at  least  with  probability  1/2  with  large  N  for  all  x  G  (Spielman  &  Srivastava,  2011). 
For  the  example  problem  considered  later,  these  bounds  are  much  tighter.  This  tells  us  that 
Lg^  holds  eigenvalues  similar  to  those  of  Lg.  These  spectrally  similar  sparse  graphs  are 
found  using  the  spectral  sparsihcation  algorithm  based  on  sampling  by  ehective  resistance 
discussed  below. 

Before  we  discuss  how  a  graph  can  be  sparsihed,  let  us  hrst  follow  the  works  of  Bollobas 
(1998)  and  Srivastava  (2010)  to  establish  an  analogy  of  a  graph  to  an  electrical  circuit.  If 
the  entire  graph  is  viewed  as  a  resistive  circuit,  we  can  dehne  a  resistance  on  the  individual 
edges  e  =  (i,  j)  of  the  graph.  According  to  Thomson’s  principle,  the  potentials  and  currents 
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in  a  resistive  circuit  distribute  themselves  so  as  to  minimize  the  total  energy  in  the  net¬ 
work  (Bollobas,  1998).  This  energy  minimization  principle  leads  to  the  concept  of  effective 
resistance. 

Effective  resistance  between  vertices  i  and  j  is  the  potential  difference  induced  between 
them  when  a  unit  current  is  injected  at  one  vertex  and  extracted  at  the  other  (Bollobas, 
1998;  Srivastava,  2010).  The  effective  graph  resistance  (also  called  as  resistance  distance) 
is  the  sum  of  the  effective  resistance  over  all  the  pairs  of  vertices  in  the  graph  Q  (Klein  & 
Randic,  1993;  Ellens  et  al.,  2011).  Rayleigh’s  monotonicity  law  states  that  pairwise  effective 
resistance  is  a  non-increasing  function  of  the  edge  weights  (Mieghem,  2011). 

In  order  to  obtain  an  expression  of  effective  resistance  for  graph  sparsihcation,  we  orient 
the  edges  of  the  original  weighted  undirected  graph  Q  with  N  vertices  and  M  edges.  We  can 
represent  any  directed  graph  by  a  signed  edge  (e)-vertex  (u)  incidence  matrix  Bg  G 
given  by 

{1  if  e  G  E  and  v  is  the  head  of  e 

—  1  if  e  G  E  and  v  is  the  tail  of  e  (6.9) 

0  otherwise. 

The  row  of  Bg  corresponding  to  an  edge  e  =  {i,j)  is  given  by  (p^  —  q^),  where  and  qj 

are  elementary  unit  vectors  in  the  i  and  j  directions,  respectively.  If  the  edge  weights  of 

the  graph  are  represented  in  a  diagonal  matrix  given  by  Cg  G  we  can  express  the 

Laplacian  matrix  based  on  the  incidence  matrix  as 

Lg  =  B^CgBg  =  ^  Wij{Pi  -  qj){p,  -  q^f.  (6.10) 

i,jeE 

The  above  relation  holds  true  only  for  undirected  graphs.  The  incidence  matrix  Bg  and  the 
diagonal  matrix  Cg  for  the  example  graph  shown  in  Fig.  6.1  (left)  are  given  by 
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For  undirected  graphs,  each  edge  is  counted  only  once  in  the  incidence  matrix  and  the 
convention  for  head  and  tail  of  an  edge  can  be  hxed  arbitrarily.  Here,  for  edge  e  =  (1,  2) 
of  the  example  graph  shown  in  Fig.  6.1  (left),  the  head  and  tail  are  considered  to  be  vertex 
1  and  2,  respectively.  The  corresponding  elementary  unit  vectors  for  the  edge  are  given  by 
Pi  =  (1,0, 0,0,0)  and  q2  =  (0,1,  0,0,0),  which  is  apparent  from  the  first  row  of  Bg  being 
Pi  -  (l2- 

For  an  edge  corresponding  to  e  =  (f,i),  the  unit  current  is  injected  at  vertex  i  and 
extracted  at  vertex  j.  Thus,  we  set  the  electrical  current  across  the  edge  to  be  (p^  —  g^). 
The  potential  induced  by  this  current  at  the  vertices  is  given  by  Lg{p^  —  qj),  where  Lg 
is  the  Moore-Penrose  pseudoinverse  of  the  Laplacian  matrix  Lg  (Srivastava,  2010).  The 
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potential  difference  across  edge  e  =  (i,  j)  is  then  given  by  (pj  —  qj)'^Lg{p^  —  qj).  Thus,  for 
unit  current,  the  effective  resistance  across  edge  e  =  {i,j)  which  corresponds  to  the  potential 
difference  can  be  expressed  as 

[Relij  =  {Pi  -  qjfL-^iPi  -  qj).  (6.11) 


The  above  expression  is  used  for  computing  effective  resistance  of  the  edges  of  the  graph.  The 
sparsihcation  of  the  original  graph  Q  =  {V,E,w}  is  performed  with  an  algorithm  Sparsify 
(Spielman  &  Srivastava,  2011)  to  produce  a  sparse  graph  Qs  =  {V,E,w}  where  w  are  the 
weights  corresponding  to  the  sparse  graph  Qs-  This  algorithm  is  based  on  the  concept  of 
effective  resistance  and  yields  a  (1  +  e)  sparse  graph  Qs-  This  sparse  graph  contains  a  reduced 
number  of  C>(iV  log(A^)/e^)  edges.  The  procedure  for  sparsihcation  Qs  =  Sparsify(^)  is 
summarized  below. 

We  hrst  create  a  list  of  the  edges  E  with  the  associated  weights  of  the  original  graph  Q. 
The  adjacency  and  Laplacian  matrices  of  graph  Q  are  constructed  from  equations  (6.1)  and 
(6.2).  The  Moore-Penrose  pseudoinverse  of  the  Laplacian  matrix  Lg  is  then  computed.  For 
each  edge  in  the  edge  list,  the  elementary  unit  vectors,  and  q^  and  edge  weights  We  are 
obtained.  Next,  the  effective  resistance  [Re\ij  corresponding  to  each  edge  is  computed  from 
equation  (6.11). 

A  random  edge  e  =  {i,j)  from  the  edge  list  of  graph  Q  is  chosen  with  probability  pe 
proportional  to  We[Re\ij-  The  edge  e  is  added  to  the  sparse  graph  Qs  with  the  weight 
given  by  We  =  We/qPe,  where  q  =  8iVlog2(iV)/e^.  We  take  integer(g)  number  of  samples 
independently  without  replacement  and  sum  the  weights  if  an  edge  is  chosen  more  than 
once.  The  resulting  graph  becomes  the  sparsihed  graph  Qs  that  satishes  equation  (6.8)  with 
eigenvalues  similar  to  those  of  Q. 

After  the  random  sampling  procedure,  a  sparsihed  adjacency  matrix  Ag^  is  obtained. 
For  convenience,  we  dehne  the  ratio  of  the  sparsihed  and  original  adjacency  matrix  weights 


as  Wij, 


’  0  otherwise. 


(6.12) 


This  sparsihcation  factor  Wij  is  related  to  the  probability  of  an  edge  e  =  {i,j)  of  graph  Q 
being  sampled.  The  edges  of  the  original  graph  that  are  not  sampled  during  the  random 
sampling  procedure  have  zero  weights  in  the  sparsihed  graph.  These  edges  are  cut  during 
sparsihcation.  The  spectral  sparsihcation  procedure  also  redistributes  the  weights  of  the 
cut  edges  among  the  other  edges  of  the  sparsihed  graph.  Thus,  cutting  of  graph  edges  is 
compensated  by  redistribution  of  the  weights  to  preserve  spectral  properties  of  the  original 
graph. 

The  spectral  sparsihcation  algorithm  described  above  produces  a  (lie)  expander  graph, 
i.e.,  sparsiher  with  strong  connectivity  properties  when  e  >  1/\/N,  where  N  is  the  number 
of  vertices  (Spielman  &  Srivastava,  2011).  With  the  concept  of  graphs  and  spectral  spar¬ 
sihcation  now  discussed,  we  consider  the  application  of  network  analysis  to  discrete  vortex 
dynamics  in  the  next  chapter. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Chapter  7 

Vortical  Interaction  Networks 


We  apply  network  analysis  to  study  the  dynamics  of  a  collection  of  point  vortices  (Saffman, 
1992;  Newton,  2001).  In  the  present  work,  the  flow  is  assumed  to  be  two-dimensional, 
incompressible,  and  inviscid.  The  spatial  domain  P  is  taken  to  be  inhnitely  large  and  the 
position  of  the  vortices  in  the  flow  field  is  denoted  by  r  =  {x,  y)  G  P.  We  consider  a  collection 
of  N  vortices  to  be  in  the  domain  V  resulting  in  the  vorticity  field  of 

N 

where  Kj  and  Vj  are  the  strength  (circulation)  and  position  of  the  j-th  vortex,  respectively, 
and  (5(-)  is  the  Dirac  delta  function.  The  motion  of  the  vortices  is  described  by  the  Biot- 
Savart  law 

dvi  Kjkx  {vj  -  Vj) 

dt  ^  2tt  \ri  —  ’ 

j=i  I  *  II 

where  k  is  the  out-of-plane  unit  normal  vector.  To  determine  the  trajectories  of  the  vortices, 
we  numerically  integrate  the  above  equation.  This  equation  is  the  basis  of  the  discrete 
vortex  methods  that  are  used  to  simulate  unsteady  vortical  flows  in  place  of  discretizing  the 
two-dimensional  Euler  equations. 

The  distance  between  the  vortices  is  dependent  on  the  location  of  the  two  point  vortices 
|rj  —  Vjl  and  the  characteristic  strength  can  be  established  by  considering  their  geometric 
mean  or  algebraic  mean  of  circulation.  We  use  the  geometric  mean  of  induced  velocity  as 
adjacency  matrix  for  the  graph  representation  of  a  set  of  point  vortices.  The  interaction  be¬ 
tween  the  system  of  discrete  point  vortices  represented  on  a  complete  graph  can  be  sparsified 
for  deriving  the  sparse  vortex  interaction  model. 

In  addition  to  considering  the  geometric  mean  of  the  circulation  for  the  weights,  algebraic 
mean  of  the  circulation  was  also  considered.  The  computation  using  the  algebraic  mean  for 
the  weights  perform  similar  to  that  with  the  above  geometric  mean.  While  we  show  that 
geometric  mean  sufficiently  captures  the  vortex  interaction,  there  may  be  some  room  for 
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Figure  7.1:  The  spatial  arrangement  of  iV  =  100  point  vortices  with  five  clusters. 


optimizing  the  choice  of  weights.  In  what  follows  we  utilize  the  geometric  mean  in  equation 
(6.7)  for  the  network  weight. 

7.1  Sparsification  of  vortex  interactions 

In  this  section,  we  construct  sparsified  representation  of  vortex  interactions  by  sparsifying 
and  redistributing  the  weights  of  the  connections  between  the  point  vortices  as  discussed  in 
§6.2.  We  consider  a  conhguration  of  N  point  vortices  in  ric  clusters.  This  setup  allows  us  to 
examine  the  effect  of  sparsification  on  a  cluster  of  vortices.  The  edge  weights  represented  in 
the  adjacency  matrix  are  given  by  equation  (6.7). 

Let  us  first  consider  a  collection  oi  N  =  100  point  vortices  with  =  5  with  each 
cluster  containing  20  vortices.  The  clusters  in  this  setup  can  be  clearly  identified  by  the 
use  of  graph  clustering  algorithms.  An  algorithm  maximizing  modularity  (Newman,  2004) 
can  be  utilized  for  the  identification  of  clusters.  One  could  also  use  /c-means  for  cluster 
identification  similar  to  the  work  performed  by  Kaiser  et  al.  (2014a).  The  point  vortices  in  a 
particular  cluster  are  given  a  normal  distribution  in  space  with  a  radial  standard  deviation  of 
(Jr  =  0.2.  The  strength  of  the  point  vortices,  within  the  individual  clusters  have  a  normal 
distribution  with  the  mean  of  k  =  0.1  and  a  standard  deviation  of  (j^  =  0.01.  The  setup 
for  this  configuration  is  shown  in  figure  7.1.  The  positions  of  the  clusters  are  normalized  by 
the  average  radial  distance  of  the  centroid  of  the  clusters  from  the  geometric  center  of  the 
system  (Ro)- 

The  original  point  vortex  distribution  modeled  as  a  complete  graph  is  shown  in  figure 

7.2  (top  left).  The  colors  of  the  point  vortices  (network  vertices)  indicate  their  unweighted 
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degree  (number  of  edge  connections).  As  the  original  graph  here  is  a  complete  graph  /Cat, 
each  vortex  is  connected  to  every  other  vortex.  Thus,  the  unweighted  degree  of  each  vortex 
in  this  example  is  A^  —  1  =  99.  The  sparsity  pattern  of  the  adjacency  matrix  for  the  complete 
graph  Aycjv  is  shown  in  hgure  7.2  (top  right).  We  also  note  in  the  hgure  the  sparsity  index 
S,  dehned  as  S'  =  n^on-zem/ with  nnon-zero  being  the  number  of  non-zero  weights  in  the 
adjacency  matrix  and  N  being  the  number  of  vortices.  As  the  diagonal  elements  of  the 
adjacency  matrix  have  zero  weights,  nnon-zero  =  A^(A^  —  1)  for  the  complete  graph.  For 
N  =  100,  a  sparsity  index  of  S  =  0.99  is  obtained  for  the  original  graph. 

We  now  perform  spectral  sparsihcation  on  the  original  graph  {Q  =  1C n)  using  the  algo¬ 
rithm  Qs  =  Sparsify(^)  described  in  §6.2  to  obtain  the  spectrally  similar  sparse  represen¬ 
tations  of  the  complete  graph  based  on  approximation  order  e.  The  corresponding  adjacency 
matrix,  Ag^  can  also  be  found.  The  vortex  network  for  sparse  graph  and  the  sparsity  pattern 
of  Agg  corresponding  to  approximation  order  of  e  =  0.5  are  shown  in  Fig.  7.2  (middle).  We 
observe  that  the  number  of  edges  (connections)  between  the  vortices  in  cluster  3  and  5  are 
reduced.  Because  the  weights  of  the  adjacency  matrix  are  inversely  proportional  to  the  dis¬ 
tance  between  the  vortices,  a  large  number  of  edges  between  the  vortices  of  the  clusters  with 
larger  distances  are  cut  during  sparsihcation.  The  sparsity  index  decreases  from  S  =  0.99 
for  the  original  complete  graph  to  S'  =  0.741  for  the  sparse  graph  (e  =  0.5),  reducing  the 
number  of  connections  by  approximately  25  percent. 

The  vortex  network  and  sparsity  patterns  for  approximation  order  of  e  =  1  are  shown  in 
hgure  7.2  (bottom).  We  observe  dense  representation  of  the  connections  between  the  clusters 
1  and  2.  Thus,  the  proximity  in  clusters  is  identihed  both  in  the  network  structure  and  the 
sparsity  patterns.  Also,  we  realize  from  the  sparsity  patterns  that  the  connections  between 
the  vortices  in  a  particular  cluster  are  maintained  while  majority  of  the  inter-cluster  ties  are 
cut  resulting  in  a  sparse  graph.  The  number  of  connections  in  the  sparse  graph  for  e  =  1  are 
reduced  dramatically  by  nearly  60  percent  with  S  =  0.4006. 

For  weighted  graphs,  the  ehective  resistance  is  computed  in  such  as  way  so  as  to  maintain 
spectral  similarity.  The  eigenspectra  {a  and  A)  of  the  adjacency  and  Laplacian  matrices  for 
the  sparse  and  original  graphs  are  compared  in  hgure  7.3.  We  observe  from  hgure  7.3  that 
the  spectra  of  the  sparse  and  original  conhgurations  are  in  good  agreement  as  expected.  The 
second  smallest  eigenvalue  of  the  Laplacian  matrix  represents  the  spectral  gap  or  algebraic 
connectivity  of  the  graph,  i.e.,  vertices  connected  by  at  least  a  single  path  (Newman,  2010). 

As  observed  from  Fig.  7.3,  the  spectral  gap  is  greater  than  zero  for  original  and  sparse 
graphs  indicating  that  the  graph  is  connected.  Thus,  sparsihcation  preserves  the  connectivity 
properties  of  the  original  graph.  The  sparsihed  graphs  slightly  underpredict  the  maximum 
eigenvalues  of  the  adjacency  and  Laplacian  matrices  compared  to  the  original  graph.  The 
preservation  of  spectral  properties  is  attractive  if  the  sparsihed  network  is  to  be  used  to 
describe  the  dynamics  of  the  vortices.  We  expect  the  overall  motion  of  the  vortices  to 
remain  similar  with  sparsihed  representations  of  the  vortex  interactions. 

Sparsihcation  helps  identify  the  individual  clusters  by  creating  subgraphs  within  a  graph 
where  the  density  of  the  edges  between  vertices  is  much  greater  than  that  outside  it.  The 
spectral  sparsihcation  procedure  leads  to  the  reduction  in  the  number  of  edges  to  0{N  log(A^)/e^) 
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Figure  7.2:  Original  and  sparsified  graphs  and  their  respective  sparsity  patterns  of  the  adja¬ 
cency  matrix  for  N  =  100  vortices.  The  color  of  the  nodes  of  the  vortex  network  represents 
the  unweighted  connections  (number  of  connected  edges)  of  the  individual  point  vortices. 
The  color  in  the  sparsity  patterns  of  the  adjacency  matrix  indicates  the  adjacency  weights 
Wij  for  original  graph  Ag  and  Wij  for  sparse  graphs  Ag^.  The  empty  white  spaces  indicate 
sparsity  in  the  adjacency  matrix.  Sparsity  index  is  S'  =  nnon-zero/-^^,  where  nnon-zero  is  the 
number  of  non-zero  weights  (elements)  in  adjacency  matrix  and  N  is  the  number  of  vortices. 
The  colored  circles  in  the  vortex  network  represent  the  individual  cluster  groups. 
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(a) 


(b) 


Figure  7.3:  Eigenvalue  spectra  of  (a)  adjacency  and  (b)  Laplacian  matrices  of  the  sparsified 
graphs  with  e  =  0.5  and  1  in  comparison  to  the  original  graph. 


for  large  N.  We  examine  the  performance  of  spectral  sparsihcation  on  the  current  example 
for  increasing  the  total  number  of  vortices.  As  in  the  example  that  was  previously  consid¬ 
ered,  we  maintain  a  constant  number  of  clusters  ric  =  5  but  increase  the  number  of  vortices 
in  each  cluster  N/ric-  The  sparsity  patterns  of  the  adjacency  matrix  for  N  =  500  and  2000 
with  approximation  order  e  =  1  are  shown  in  Fig.  7.4.  We  observe  that  as  the  number  of 
vortices  increases,  the  adjacency  matrix  of  the  sparsihed  graph  becomes  increasingly  sparse. 
Similar  to  the  case  with  N  =  100  vortices,  the  majority  of  the  ties  within  a  cluster  are  main¬ 
tained  while  a  large  number  of  inter-cluster  ties  are  cut.  The  sparsity  index  decreases  from 
S  =  0.1784  to  0.0707  for  N  =  500  to  2000,  which  is  a  substantial  amount  of  sparsihcation. 
Let  us  further  show  the  sparsity  index  S  for  increasing  N  with  different  approximation  order 
e  in  Fig.  7.5.  The  sparsity  index  S  decreases  with  an  increase  in  e  and  the  number  of  vor¬ 
tices.  The  expected  behavior  of  the  sparsity  index  S  =  0{N  \og{N)) / [N"^)  =  0(\og{N)/N) 
is  observed  for  larger  N .  The  trend  deviates  from  the  expected  behavior  for  lower  N  as  the 
availability  of  the  connections  for  redistribution  of  the  weights  of  the  sparsihed  connections 
decreases. 

We  observe  that  the  sparsihcation  algorithm  provides  us  with  a  heavily  sparsihed  model  to 
produce  a  computationally  tractable  representation  of  the  vortex-to-vortex  interaction.  We 
further  note  that  the  sparsihcation  algorithm  can  be  easily  parallelized,  which  is  attractive 
to  further  reduce  the  computational  wait  time.  The  bulk  motion  of  clusters  of  point  vortices 
could  be  tracked  ehectively  by  tracking  their  individual  centroids.  The  present  approach 
may  lead  to  algorithms  similar  to  fast  particle  summation  methods  where  near  and  far-held 
ehects  are  taken  into  consideration  (Greengard  &  Rokhlin,  1987). 
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Figure  7.4:  Sparsity  patterns  of  adjacency  matrix  for  N  =  500  and  2000  for  approximation 
order  e  =  1.  The  color  in  the  sparsity  patterns  of  the  adjacency  matrix  denotes  the  adjacency 
weights  Wij  with  the  empty  white  space  indicating  sparsity. 
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Figure  7.5:  Variation  of  sparsity  index  S  for  sparsihcation  with  increasing  number  of  point 
vortices.  The  reference  black  line  represents  log(V)/iV. 
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7.2  Sparsified-dynamics  model 

The  sparsified  vortex  interactions  can  be  used  for  analyzing  the  dynamics  of  the  system  of 
discrete  point  vortices.  Spectral  sparsihcation  can  be  viewed  as  adjusting  the  strength  of  the 
vortices  in  the  Biot-Savart  law  through  the  sparsihcation  factor  W  given  by  equation  (6.12). 
We  can  incorporate  the  weights  of  the  sparse  graph  into  the  Biot-Savart  law  as 


(7,3) 


where  Wij  is  the  sparsihcation  factor.  The  above  expression  signihcantly  reduces  the  amount 
of  computation  compared  to  the  original  Biot-Savart  law,  because  the  sparsihcation  has  made 
a  signihcant  number  of  elements  of  W  to  be  zero,  as  seen  in  hgures  7.2  and  7.4. 

For  evaluating  the  ehect  of  sparsihcation  on  dynamics  of  the  point  vortices,  we  consider 
the  same  example  presented  in  §7.1  with  N  =  100  vortices  as  the  initial  condition.  The 
objective  here  is  to  accurately  capture  the  dynamics  of  the  vortex  clusters  with  the  sparsihed 
Biot-Savart  law  given  by  equation  (7.3).  The  bulk  motions  of  the  clusters  of  point  vortices 
are  well  represented  by  the  centroid  of  the  individual  clusters.  The  original  and  the  sparsihed 
dynamics  are  given  by  equations  (7.2)  and  (7.3),  respectively,  and  are  integrated  in  time  with 
the  fourth-order  Runge-Kutta  method.  The  results  from  the  sparsihed-dynamics  model  are 
shown  in  hgures  7.6  and  7.7. 

For  the  present  analysis,  time  is  non-dimensionalized  by  considering  the  cumulative  cir¬ 
culation  of  the  vortices  of  the  cluster  (F)  and  average  radial  distance  of  the  centroid  of  the 
clusters  from  the  geometric  center  of  the  overall  system  at  the  initial  time  {Rq)-  As  the 
advective  velocity  of  an  isolated  cluster  of  point  vortices  of  strength  F  at  Ro  is  given  by 
u*  =  F/27r i?o,  the  non-dimensional  time  can  be  deduced  as  tu* / Ro  =  tV /2'kR^^.  In  the  cur¬ 
rent  example,  each  individual  cluster  has  a  collection  of  vortices  with  their  strength  having 
a  mean  of  /t  =  0.1  (with  a  normal  distribution).  Thus,  for  n  vortices  in  a  cluster,  F  =  nR. 
The  error  in  position  of  the  centroid  of  the  clusters  for  the  sparsihed  setting  with  respect  to 
the  original  setting  given  by  \r^  —  r\  is  non-dimensionalized  by  the  average  radial  distance  of 
the  centroid  of  the  clusters  from  the  geometric  center  of  the  system  at  the  initial  time  {Ro)- 

Let  us  consider  the  case  where  sparsihcation  is  performed  only  once  to  determine  the 
sparsihcation  factor  Wij  before  initiating  the  time  integration.  The  same  sparsihcation  factor 
is  used  throughout  the  time  integration.  This  procedure  is  inexpensive  as  sparsihcation  is 
performed  only  once  a  priori.  The  trace  of  the  positions  of  the  centroid  of  the  individual 
clusters  and  their  error  for  approximation  orders  of  e  =  0.5  and  1  for  a  single  sparsihcation  are 
shown  in  hgures  7.6  (top)  and  7.7  (top),  respectively.  Considering  the  number  of  connections 
cut  between  the  vortices  for  sparse  approximations,  the  trajectories  of  the  centroid  of  the 
clusters  based  on  sparse  vortex  dynamics  agrees  reasonably  with  those  from  full  dynamics. 
On  comparing  the  errors,  we  hnd  that  the  errors  with  approximation  order  of  e  =  0.5  is  less 
than  those  with  e  =  1,  which  is  expected.  Despite  the  Biot-Savart  law  being  nonlinear,  the 
present  sparsihed  approach  achieves  a  reasonable  agreement  with  the  full  nonlinear  solution. 
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If  the  clusters  of  vortices  are  far  from  each  other,  the  connections  cut  between  their 
vortices  do  not  cause  a  large  change  in  the  dynamics  while  the  connections  cut  between 
clusters  close  to  each  other  affect  the  dynamics  considerably.  As  time  progresses,  the  relative 
distance  between  vortices  change,  resulting  in  an  increase  in  error.  For  tV /2'kR^  >  1,  we 
observe  that  the  errors  in  sparse  vortex  dynamics  increases  to  ~  0.5i?o.  As  the  adjacency 
weights  Wij  and  the  sparsification  factor  Wij  are  based  on  the  initial  position  of  the  vortices, 
the  ratio  does  not  hold  over  large  times.  This  calls  for  adjustment  of  the  weights  and  periodic 
sparsihcation  to  adapt  to  the  dynamically  changing  position  of  the  vortices. 

We  thus  consider  resparsihcation  of  the  graph  representing  the  system  of  point  vortices 
periodically  at  {Q.ltV /2ttR^^)  and  (0.01tr/27ri?Q)  for  both  approximation  orders.  The  trace 
of  the  centroids  and  error  for  approximation  orders  of  e  =  0.5  and  1  with  resparsihcation  are 
shown  in  middle  and  bottom  subhgures,  respectively,  of  hgures  7.6  and  7.7.  We  observe  that 
the  error  decreases  signihcantly  with  resparsihcation.  The  error  in  position  reduces  from 
0(10“^)  to  0(10“^)  with  resparsihcation.  We  notice  that  the  centroid  trajectories  based 
on  the  sparse  and  original  calculations  become  increasingly  similar  with  resparsihcation. 
Resparsihcation  updates  the  sparsihcation  factor  periodically  based  on  the  position  and 
strength  of  the  vortices  and  decreases  the  error  by  an  order  of  magnitude.  Thus,  the  nonlinear 
evolution  of  discrete  vortex  dynamics  is  well  predicted  by  the  sparsihcation  techniques. 

One  of  the  advantages  of  sparsihcation  is  the  decreased  computational  cost  due  to  in¬ 
creased  sparsity  of  the  connections  between  the  vortices.  This  could  potentially  lead  to  design 
of  faster  algorithms  based  on  sparsihcation  strategies.  Let  us  hrst  evaluate  the  offline  cost  of 
computing  ehective  resistance  and  random  sampling  required  for  spectral  sparsihcation.  The 
time  required  for  computing  ehective  resistance  R  and  time  required  for  random  sampling 
tg  is  shown  in  hgures  7.8(a)  and  (b)  respectively.  All  computations  were  performed  in  MAT- 
LAB  on  an  iMac  with  3.4  GHz  Intel  Core  i7  processor.  We  can  observe  that  for  each  edge, 
the  computation  time  of  ehective  resistance  is  (!7(log(A^)),  hence  requiring  0{N‘^  \og{N))  for 
the  overall  ehective  resistance  computation.  The  random  sampling  procedure  takes  0{N‘^) 
time  for  all  the  approximation  orders.  We  compare  the  computational  time  required  at 
each  step  of  numerical  integration  of  the  Biot-Savart  law  for  the  original  and  sparse  conhg- 
uration.  We  can  observe  from  hgure  7.8(c)  that  the  original  conhguration  takes  0{N‘^)  time 
equivalent  to  the  number  of  edges  in  the  complete  graph,  while  for  the  sparse  conhguration 
the  computation  time  is  reduced  to  C>(A^ log(A^)). 

There  is  always  a  tradeoh  amongst  the  level  of  sparsihcation,  the  error  that  appears  in 
the  dynamics,  and  the  associated  computational  cost.  The  sparser  the  network  is,  the  faster 
the  computation  can  become  but  may  compromise  accuracy.  To  increase  the  computational 
accuracy,  resparsihcation  can  be  performed  but  can  introduce  an  additional  computational 
load.  We  evaluate  the  computational  cost  for  numerical  integration  until  a  total  time  of 
0.1fr/27ri?Q.  The  time  required  for  resparsihcation  and  numerical  integration  for  the  two 
diherent  resparsihcation  frequencies  (performed  every  0.1tr/27iRl  and  Q.QltV /2ttR^)  con¬ 
sidered  in  this  work  for  original  and  sparse  conhguration  with  e  =  1  is  shown  in  hgure  7.9. 
It  should  be  noted  that  for  the  resparsihcation  performed  every  0.1fr/27ri?Q,  resparsihca¬ 
tion  is  performed  once  for  the  total  time  considered.  For  resparsihcation  conducted  every 
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Figure  7.6:  Trajectories  of  the  centroids  (left)  and  error  in  position  of  centroids  (right) 
of  the  vortex  clusters  based  on  sparse  graph  (e  =  0.5)  and  original  graph  with  time  for 
N  =  100,  K  =  0.1.  (a)  trajectories  and  error  based  on  sparsihcation  performed  only  at  the 
initial  step  while  (b)  results  based  on  resparsihcation  performed  at  every  (0.1tr/27ri?o)  and 
(c)  results  based  on  resparsihcation  performed  at  every  /2ttR^).  The  full  nonlinear 

solution  is  shown  with  dashed  lines. 
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Figure  7.7:  Trajectories  of  the  centroids  (left)  and  error  in  position  of  centroids  (right) 
of  the  vortex  clusters  based  on  sparse  graph  (e  =  1)  and  original  graph  with  time  for 
N  =  100,  K  =  0.1.  (a)  trajectories  and  error  based  on  sparsihcation  performed  only  at  the 
initial  step  while  (b)  results  based  on  resparsihcation  performed  at  every  (0.1tr/27ri?o)  and 
(c)  results  based  on  resparsihcation  performed  at  every  /2ttR^).  The  full  nonlinear 

solution  is  shown  with  dashed  lines. 
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(b) 


Figure  7.8:  The  time  required  for  computation  of  (a)  effective  resistance  and  (b)  random 
sampling  for  different  approximation  orders  for  one-time  sparsification.  The  dashed  lines 
indicate  the  expected  trends,  (c)  Time  required  for  numerical  integration  at  each  step  for 
original  and  sparsihed  conhguration  with  e  =  1. 
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Figure  7.9:  The  time  required  for  numerical  integration  for  a  total  time  of  O.lfF /2ttR^  with 
different  resparsihcation  frequencies. 


0.01fr/27ri7o5  resparshcation  is  performed  ten  times.  We  can  see  that  the  time  required  for 
resparsihcation  at  every  O.ltF is  considerably  less  compared  to  the  original  conhgura- 
tion.  With  increase  in  resparsihcation  frequency,  though  the  time  required  is  less  compared 
to  the  original  conhguration,  there  is  increased  cost  associated  with  resparsihcation. 

To  illustrate  that  the  trajectory  prediction  using  sparsihcation  cannot  be  performed 
naively,  we  compare  the  results  from  spectral  sparsihcation  with  those  based  on  random 
edge  removal  (without  redistribution  of  edge  weights).  We  keeep  the  number  of  cuts  identi¬ 
cal  to  that  achieved  by  spectral  sparsihcation  for  fairness.  The  sparsity  patterns  for  random 
cuts  with  sparsity  S  =  0.7410  and  0.4006  and  the  corresponding  Laplacian  eigenspectra 
are  shown  in  hgure  7.10.  We  observe  that  the  spectra  from  the  randomly  sparsihed  graphs 
grossly  under-predicts  that  of  the  original  graph. 

The  trace  of  the  positions  of  the  centroid  of  the  individual  clusters  and  their  error  for 
random  cuts  with  sparsity  S  =  0.7410  and  0.4006  for  one-time  sparsihcation  are  shown 
in  hgure  7.11.  The  random-cut  approximation  performs  poorly  as  compared  to  spectral 
sparsihcation  and  the  dynamics  of  the  centroid  of  the  vortex  clusters  are  not  captured.  As 
the  connections  between  the  vortices  are  randomly  cut,  the  centroids  of  the  clusters  move 
slower  compared  to  the  original  conhguration  as  expected.  This  is  due  to  the  loss  of  induced 
velocity  with  the  absence  of  weight  redistribution  from  random  sparsihcation.  Spectral 
sparsihcation,  on  the  other  hand,  redistributes  the  weights  to  prevent  the  loss  of  the  overall 
interactions  amongst  the  set  of  point  vortices. 

It  is  well  known  that  fast  multipole  methods  (Greengard  &  Rokhlin,  1987)  can  also  reduce 
the  0{N‘^)  velocity  evaluation  signihcantly.  The  particle-box  and  box-box  schemes  reduce 
the  computational  complexity  from  0{N‘^)  to  0{N  \og{N))  and  0{N),  respectively.  As  seen 
in  this  section,  spectral  sparsihcation  reduces  the  computational  load  to  0{N  log{N))  for 
computing  the  velocity  of  the  vortices.  Fast  multipole  methods  approximate  the  ehect  of 
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Figure  7.10:  Sparsity  patterns  of  the  adjacency  matrix  for  N  =  100  vortices  with  random 
cuts  with  sparsity  (a)  S  =  0.7410  and  (b)  S  =  0.4006.  The  color  in  the  sparsity  patterns 
of  the  adjacency  matrix  indicates  the  adjacency  weights  Wij  for  S  =  0.7410  and  0.4006.  (c) 
Eigenvalue  spectra  of  Laplacian  matrices  of  the  randomly  sparsihed  graphs  with  S  =  0.7410 
and  0.4006  in  comparison  to  the  original  graph.  The  legend  for  the  sparsity  patterns  is 
similar  to  that  of  sparsity  patterns  in  hgure  7.2. 
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Figure  7.11:  (a)  Trajectories  of  the  centroids  and  (b)  error  in  position  of  centroids  (right) 
of  the  vortex  clusters  based  on  random  cuts  with  sparsity  S  =  0.7410  and  0.4006  for  N  = 
100,  it  =  0.1. 


cluster  of  particles  at  a  certain  distance  by  multipole  expansions  and  organizing  the  parti¬ 
cles  to  a  hierarchy  of  clusters  (Cottet  &  Koumoutsakos,  2000).  In  order  to  compensate  for 
the  sparsihed  connections  in  spectral  sparsihcation,  the  weights  are  redistributed  among  the 
other  edge  connections  to  compute  the  sparsihed  dynamics.  In  an  analogy  to  fast  multipole 
methods,  the  interaction  list  depends  on  the  connections  that  have  direct  impact  on  the  vor¬ 
tices  and  computational  accuracy  may  degrade  when  the  vortices  move  in  space.  Similar  to 
reconstruction  of  tree  data  structures  in  fast  multipole  methods,  performing  resparsihcation 
periodically  increases  the  accuracy  by  re-evaluating  the  associated  weights  as  the  positions 
of  the  vortices  evolve  over  time.  We  note  that  even  with  sparsihcation,  invariants  of  discrete 
vortices  are  conserved  as  discussed  below. 
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7.2.1  Conservation  properties 

For  a  discrete  set  of  point  vortices,  angular  impulse  and  linear  impulse  are  among  the  con¬ 
served  quantities.  The  coordinates  of  center  of  vorticity  obtained  from  linear  impulse  and 
the  length  of  dispersion  about  the  center  of  vorticity  obtained  from  angular  impulse,  when 
the  strength  of  vortices  is  of  the  same  sign,  remain  constant.  The  center  of  vorticity  {X,Y) 
and  length  of  dispersion  of  vorticity  D  are  given  by 


2^j=l 


L^i=l 

2^j=l 


Y  = 


K  ■  '' 

z^i=i 


(7.4) 

(7.5) 


where  Ki  is  the  strength  of  the  i-th  vortex  and  {x,  y)  are  the  vortex  positions  (Batchelor,  2000; 
Newton,  2001).  The  circulation  for  n  vortices  in  a  cluster  is  T  =  nn.  The  total  circulation  of 
the  system  of  Uc  =  5  clusters  is  given  by  T^  =  =  ’^cT.  Another  conserved  quantity 

is  the  Hamiltonian.  The  Hamiltonian  is  dehned  as  the  interaction  energy  of  the  system  of 
point  vortices  and  is  given  by 


^  N 

^  =  log  \/ (Xi  -  +  iVi  -  VjY-  (7.6) 

We  non-dimensionalize  the  Hamiltonian  using  the  total  circulation  of  ric  clusters  and  the 
average  radial  distance  of  centroid  of  the  clusters  at  initial  time  Ro- 

We  compare  the  error  in  the  Hamiltonian  {H^{r/Ro)—H{r/Ro))/T^,  the  center  of  vorticity 
{{X^  —  X)/ Ro,  {Y^  —  Y)/Ro),  and  the  square  of  length  of  dispersion  of  vorticity  {D‘^  —  D‘^)/Rl 
for  approximation  orders  of  e  =  0.5  and  1  over  time  in  hgure  7.12.  The  preservation  of  center 
of  vorticity  and  length  of  dispersion  implies  the  conservation  of  linear  and  angular  impulse. 
Here,  the  variables  with  subscript  e  denote  those  based  on  a  sparsihed  model.  We  observe 
that  the  error  in  the  Hamiltonian  is  of  0(10“^)  for  e  =  0.5  and  0(10“^)  for  e  =  1.  The 
error  in  the  square  of  dispersion  is  of  (T(10“®)  while  the  error  in  the  center  of  vorticity  is  of 
0(10“®)  with  sparsihcation.  While  not  shown,  resparsihcation  of  these  vortices  periodically 
performed  at  every  (O.ltT /27iRl)  and  /2'kR^)  maintain  similar  error  levels  for  the 

invariants.  We  also  note  that  the  circulation  of  each  vortices  has  not  been  altered,  which 
conserves  the  individual  and  overall  circulation  over  time. 

As  seen  in  §6.2,  sparsihcation  by  effective  resistance  is  based  on  the  energy-minimization 
principle.  The  effective  resistance  is  dehned  such  that  the  total  energy  of  the  system  re¬ 
mains  constant.  This  is  rehected  well  in  the  conservation  of  the  Hamiltonian  representing 
the  interaction  energy  of  the  discrete  point  vortices.  In  addition,  sparsihcation  preserves 
other  invariants  including  linear  and  angular  impulse  of  the  discrete  set  of  point  vortices. 
We  note  in  passing  that  even  with  resparsihcation,  the  invariants  were  observed  to  be  con¬ 
served.  The  conservation  of  the  invariants  for  the  sparse  conhgurations  further  encourages 
the  applicability  of  sparsihcation  strategies  on  discrete  vortex  dynamics. 
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Figure  7.12:  Error  in  the  (a)  Hamiltonian,  (b)  square  of  length  of  dispersion,  (c)  x-position 
of  center  of  vorticity  and  (d)  ^/-position  of  center  of  vorticity  based  on  single  sparsification 
for  N  =  100,  R  =  0.1  and  e  =  {0.5, 1}. 
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7.2.2  Comparison  with  Reduced-order  model 

We  comment  on  the  similarities  and  differences  between  the  present  sparsihed-dynamics 
model  and  the  reduced-order  model.  Both  models  share  the  general  objective  of  deriving 
a  model  that  captures  the  full-order  physics  in  a  distilled  manner.  Reduced-order  models 
achieve  such  goals  by  reducing  the  dimension  of  the  state  variable.  One  such  technique  is 
the  Galerkin  projection  of  the  Navier-Stokes  equations  using  a  set  of  spatial  bases  (Noack 
et  al,  2005;  Rowley  et  al,  2004),  such  as  those  determined  from  POD  (Berkooz  et  al, 
1993;  Holmes  et  al.,  1996).  One  can  further  consider  incorporating  the  effect  of  input  and 
output  dynamics  by  using  balanced  truncation  or  the  eigensystem  realization  algorithm 
(Rowley,  2005;  Ma  et  al,  2009).  The  resulting  reduced-order  models  are  generally  described 
by  ordinary-differential  equations  for  the  temporal  coefficients  with  reduced  dimensions. 

On  the  other  hand,  the  present  sparsihed-dynamics  model  does  not  reduce  the  dimen¬ 
sionality  of  the  state  variable.  We  instead  focus  on  reducing  the  number  of  connections  used 
to  capture  the  overall  dynamics  of  the  huid  how.  The  reduction  in  the  number  of  connections 
is  performed  based  on  the  concepts  from  network  analysis  and  graph  theory.  Although  we 
cut  a  large  number  of  the  edges  in  the  graph  representation  of  the  dynamical  interaction, 
we  redistribute  the  weights  associated  with  the  edges  to  maintain  properties  of  the  graph. 
This  procedure  is  based  upon  approximating  the  spectral  properties  of  the  graph  and  does 
not  require  selecting  the  spatial  basis  functions  unlike  the  Galerkin-projection  based  models. 
One  nice  feature  of  the  sparsihed-dynamics  model  is  its  ability  to  conserve  physical  vari¬ 
ables  such  as  the  Hamiltonian,  angular  impulse,  and  linear  impulse,  as  discussed  previously. 
Hence,  the  sparsihed  dynamics  is  able  to  predict  the  full  dynamics  as  demonstrated  by  the 
example  with  discrete  vortices.  By  highlighting  the  interactions  amongst  a  set  of  vortices, 
we  are  able  to  determine  which  interactions  amongst  the  vortices  are  important  in  guiding 
the  overall  motion  of  vortices. 

We  believe  the  sparsihed-dynamics  model  has  promising  potential  to  model  various  types 
of  huid  how  by  considering  the  modal  structures  as  a  abstraction  of  graph  nodes.  The 
application  of  sparsihed  dynamics  models  towards  how  control  problems  may  also  become 
fruitful  as  the  the  computational  time  necessary  to  capture  the  complex  behavior  of  the  how 
is  reduced  and  the  interaction  of  nodes  or  how  structures  are  well-captured,  which  has  been 
a  lacking  feature  in  linear  dynamics  model  and  linear  stability  analysis.  We  anticipate  that 
the  presently  proposed  model  can  lead  to  potential  feedback  control  of  huid  how  (Bagheri 
et  al,  2009;  Ahuja  &  Rowley,  2010)  but  with  nonlinear  interactions  emphasized.  One  of  the 
open  questions  in  extending  the  sparsihed  dynamics  model  for  a  wide  variety  of  huid  how 
problems  is  the  choice  of  the  variables  or  modes  to  be  used  for  graph  nodes.  This  issue  is 
currently  being  examined  and  will  be  reported  in  upcoming  studies.  In  the  next  chapter,  we 
extend  the  vortical  interaction  analysis  to  turbulent  hows  using  network  theory. 
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Chapter  8 

Turbulent  Interaction  Networks 


In  this  study,  we  focus  on  unforced  two-dimensional  isotropic  turbulence  in  a  periodic  box 
and  assess  the  influence  of  the  vorticity  distribution  over  a  Cartesian  domain.  Here,  the 
two-dimensional  vorticity  field  reduces  to  uj{x,t)  =  u{x,y,t)ez  with  denoting  the  unit 
normal  plane  vector.  Modeling  the  vortical  component  for  each  discrete  Cartesian  element 
as  a  line  vortex,  we  can  evaluate  how  fluid  elements  influence  each  other,  as  depicted  in 
Fig.  8.1.  We  define,  Ki  =  u}{xi)AxAii  is  the  circulation  of  fluid  element  i  with  side  lengths  of 
Ax  and  Ay.  The  superposition  of  the  induced  velocity  from  all  other  fluid  elements  provides 
the  advective  velocity  of  the  fluid  element.  Detailed  discussions  on  using  point  vortices  to 
develop  the  network-theoretic  framework  for  describing  unsteady  vortical  flows  can  be  found 
in  Nair  &  Taira  (2015)  and  previous  chapter.  We  note  in  passing  that  adjacency  matrices 
are  commonly  defined  with  positive  weights  as  considered  here,  but  they  can  be  relaxed  to 
accommodate  positive  and  negative  weights  within  the  context  of  vortical  interactions.  This 
point  will  be  revisited  later. 

In  the  present  study,  the  influence  from  the  neighboring  periodic  vortex  images  are  also 
accounted  for  in  the  analysis.  This  formulation  yields  a  full  matrix  except  for  its  diagonal 
entries  that  are  identically  zero.  In  assessing  the  strength  of  the  vortical  interaction  between 
two  fluid  elements,  we  utilize  Eq.  (6.6)  to  perform  network  analysis  to  extract  the  spatial 
connectivity  structure.  This  approach  has  been  successful  in  capturing  the  nonlinear  vortex 
dynamics  and  modeling  the  trajectories  of  vortex  clusters  (Nair  &  Taira,  2015).  The  average 
induced  velocity  is  used  as  the  adjacency  matrix.  Note  that  the  geometric  mean  can  be 
alternatively  chosen  and  yields  similar  results.  In  general,  the  adjacency  matrix  can  be 
formulated  in  an  asymmetric  manner: 

+  (1  -  (j))uj^i 

0  otherwise 

Here  the  parameter  0  takes  a  value  between  0  and  1.  For  the  aforementioned  symmetric 
formulation  in  Eq.  (6.6),  0  is  selected  as  1/2.  When  0  =  0  and  1,  the  adjacency  matrix  Aij 
are  defined  by  the  velocity  imposed  to  the  other  elements  {Aij  =  Uj^i)  and  upon  themselves 
{Aij  =  Ui^j),  respectively,  for  i  ^  j.  We  mainly  focus  on  the  use  of  the  symmetric  adjacency 
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Figure  8.1:  Interaction  of  fluid  elements  in  two-dimensional  turbulence.  The  strength  of  the 
vortical  interaction  between  elements  i  and  j  having  vorticity  cuj  and  Uj  is  quantihed  through 
the  induced  velocities  Ut^j  and  Uj^i,  respectively.  For  discretizing  the  Cartesian  domain,  we 
take  rix  and  Uy  points  in  the  horizontal  and  vertical  directions,  respectively,  providing  the 
adjacency  matrix  A  of  size  n  x  n  with  n  =  UxUy.  Shown  in  the  background  with  a  contour 
plot  is  the  corresponding  vorticity  held  with  initial  i?e(fo)  =  814  at  t  =  18. 

matrix  in  this  work  but  will  consider  the  asymmetric  formulation  briehy  to  highlight  the 
difference  from  a  physical  point  of  view  in  the  next  section.  We  note  in  passing  that  the 
theoretical  tools  for  symmetric  adjacency  matrices  are  more  widely  available  compared  to 
the  asymmetric  matrices. 

The  how  held  analyzed  in  this  study  is  obtained  from  direct  numerical  simulation  on  a 
square  bi-periodic  computational  domain  (x,  y)  G  [0,  L]  x  [0,  L]  with  a  grid  size  of  nix  x 
niy  =  1024  X  1024.  The  unforced  two-dimensional  incompressible  isotropic  turbulent  how  is 
simulated  by  numerically  solving  the  two-dimensional  vorticity  transport  equation 

du  du  _  1  d'^oj 

dt  dxj  Re  dxjdxj  ’ 

where  u  and  uj  are  the  velocity  and  vorticity  variables,  respectively.  The  simulation  is  per¬ 
formed  with  the  Fourier  spectral  method  and  the  fourth-order  Runge-Kutta  time  integration 
scheme  (Canuto  et  ai,  1988).  The  vorticity  held  is  initialized  with  a  smooth  distribution 
comprised  of  a  large  number  100)  of  superposed  vortices  (Taylor,  1918)  with  random 
strengths,  core  sizes,  and  locations.  The  initial  core  sizes  are  selected  to  be  sufficiently  small 
compared  to  the  size  of  the  computational  domain  (McWilliams,  1984)  arranged  in  random 
positions.  The  velocity  variable  is  normalized  by  the  square  root  of  the  spatial  average  of 
the  initial  kinetic  energy  u*{to)  =  [^^(to)]^'^^,  where  the  overline  denotes  the  spatial  average. 
The  spatial  length  and  time  scales  are  non-dimensionalized  by  the  initial  integral  length  scale 
I*  (to)  =  [2u‘^(to)/u)^(to)Y''‘^  and  the  initial  eddy  turnover  time  =  I*  (to) / u*  (to) ,  respectively. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


The  Reynolds  number  is  defined  accordingly  as  Re  =  u*l* /u  where  u  is  the  kinematic  vis¬ 
cosity.  In  this  study,  turbulent  fiows  with  initial  Reynolds  numbers  of  Re(to)  =  75,  439,  814, 
1607,  and  2485  are  selected. 

8.1  Network-based  characterization 

We  identify  the  underlying  network  structure  and  characteristics  of  two-dimensional  turbu¬ 
lence  based  on  the  aforementioned  symmetric  adjacency  weights.  The  time-evolving  vorticity 
field  is  obtained  from  a  two-dimensional  incompressible  bi-periodic  direct  numerical  simula¬ 
tion  (Canuto  et  al.,  1988)  for  unforced  isotropic  turbulence.  Given  the  vorticity  field  over  a 
Cartesian  grid,  each  fluid  element  is  considered  to  be  connected  to  all  other  elements  through 
vortical  network  edges.  The  resulting  fluid  flow  network  can  in  fact  be  described  by  a  com¬ 
plete  graph  with  a  range  of  weights.  Next,  we  visualize  the  network  edges  with  transparent 
gray  scale  corresponding  to  the  adjacency  weight,  as  shown  in  Fig.  8.2(a).  The  captured 
structure  reveals  the  turbulent  network.  Some  regions  in  the  flow  have  a  large  number  of 
strong  connections  corresponding  to  larger  stronger  vortices  seen  in  red,  serving  as  primary 
network  hubs.  Note  that  these  strong  vortices  induce  velocities  over  long  distances.  Mod¬ 
erate  size  vortices  that  act  as  secondary  hubs  also  possess  dominant  connections  to  primary 
hubs  and  other  secondary  hubs.  In  contrast,  fluid  elements  corresponding  to  smaller,  weaker 
eddies,  shown  in  blue,  generally  have  influence  only  in  their  vicinity.  The  node  strength  dis¬ 
tribution  (sj  =  Aij)  over  space  shows  that  the  vortices  with  large  circulation  have  larger 
strength,  as  illustrated  in  Fig.  8.2(b).  The  node  strength  distribution  over  space  enables 
us  to  distinguish  secondary  and  primary  hubs,  which  may  not  be  easily  differentiated  from 
simply  visualizing  the  vorticity  field  or  the  Q  criterion  in  a  traditional  manner.  For  instance, 
see  the  green  vortices  in  (b)  which  can  appear  similar  to  primary  ones  in  vorticity  level. 

Plotting  the  probability  of  strength  distribution  P{s)  over  the  strength  s  of  fluid  elements 
in  Fig.  8.2(c),  we  find  that  two-dimensional  isotropic  turbulence  network  has  a  power-law 
distribution  P{s)  ~  with  7  =  2.7  at  the  time  shown.  This  tells  us  that  the  vortex  inter¬ 
actions  in  turbulence  can  be  characterized  by  a  weighted  scale-free  network.  This  realization 
enables  the  interaction-based  analysis  of  turbulent  fiows  from  a  new  perspective  through  net¬ 
work  theory  (Newman,  2010;  Cohen  &  Havlin,  2010).  In  particular,  this  type  of  network  is 
known  to  have  certain  resilience  properties  as  we  will  explore  later  in  this  section.  Also  shown 
in  Fig.  8.2(c)  in  gray  are  the  degree  distributions  for  asymmetric  adjacency  formulations. 
The  out  and  in-degree  distributions  can  be  found  by  setting  0  =  0  and  1,  respectively,  in 
Eq.  (8.1).  It  can  be  observed  that  the  scale-free  symmetric  distribution  is  mostly  comprised 
of  the  out-degree  components,  which  describe  how  each  vortical  element  influences  all  other 
elements  (i.e.,  Uj^i).  In  contrast,  we  find  that  the  in-degree  distribution  has  a  single  peak 
which  conveys  that  all  fluid  elements  receive  a  similar  amount  of  collective  influence  from 
vortices  in  the  flow  field.  We  have  found  that  the  scale-free  property  of  two-dimensional 
isotropic  turbulence  is  most  well-captured  by  the  symmetric  weights  compared  to  the  other 
asymmetric  formulations.  It  is  also  possible  to  examine  the  strength  distribution  taking 
positive  and  negative  values  of  circulations,  as  we  have  briefly  discussed  in  Section  2.  Uti- 
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Figure  8.2:  The  scale-free  network  of  vortical  interactions  in  two-dimensional  turbulence 
with  initial  Re{to)  =  814.  (a)  Turbulent  network  structure  overlaid  on  the  vorticity  field 
with  the  darkness  of  the  network  edges  corresponding  to  the  values  of  the  adjacency  weights 
{t  =  18).  (b)  Contour  plot  of  the  node  strength  s  distribution.  Vortex  cores  having  high 
degree  of  connectivity  act  as  hubs  in  the  turbulent  vortical  network,  (c)  The  corresponding 
node  strength  probability  distribution  exhibiting  the  scale-free  characteristics  with  P  ~ 

The  same  contour  level  is  shared  by  (b)  and  (c).  Also  shown  in  the  background  of  (c)  in  gray 
are  the  out  and  in-degree  distributions  (0  =  0  and  1,  respectivley) .  The  network  visualized 
in  (a)  does  not  show  interactions  from  periodic  images  and  uses  32  x  32  nodes  for  graphical 
clarity. 
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Figure  8.3:  The  dynamics  of  turbulent  network  with  Re{to)  =  814.  (a)  Kinetic  energy  and 
(b)  strength  distribution  of  two-dimensional  isotropic  turbulence  for  t  =  15,  30,  75,  150,  and 
300  (line  colors  represent  time).  The  inset  plots  in  (b)  show  the  corresponding  vorticity  helds. 
The  kinetic  energy  E{k)  is  shown  over  the  wavenumber  k  exhibiting  the  asymptotic  prohle 
of  E{k)  ~  k~^.  The  strength  distribution  P{s)  displays  the  scale-free  property  of  P{s)  ~  s~'^ 
over  node  strength  s.  (c)  The  corresponding  exponents  7,  71,  and  72  are  shown.  Later  in 
time  the  strength  distribution  exhibits  the  emergence  of  two  distributions,  P{s)  ~  and 


lizing  positive  and  negative  weights,  their  strength  distribution  can  also  exhibit  a  scale-free 
behavior  but  with  network  strength  having  both  negative  and  positive  values.  This  leads 
to  a  symmetric  strength  distribution  over  the  strength  with  resemblance  to  the  probability 
density  function  of  scaled  displacements  (Weiss  et  al,  1998).  In  what  follows,  results  based 
on  the  symmetric  adjacency  matrix  (using  the  magnitude  of  induced  velocity)  are  presented. 

Let  us  further  examine  the  time- varying  properties  of  the  turbulent  network.  In  unforced 
turbulence,  the  kinetic  energy  of  the  flow  decreases  over  time  due  to  viscous  dissipation  as 
shown  in  Fig.  8.3(a).  The  strength  distribution  P{s)  of  the  turbulence  network  and  the 
corresponding  flow  held  snapshots  are  presented  in  Fig.  8.3(b).  Turbulent  how  is  comprised 
of  vortical  structures  over  a  wide  range  of  spatial  scales  initially.  The  distribution  P{s) 
exhibits  scale-free  characteristics  for  t  <  30  with  P{s)  ~  s~'^,  where  7  ~  2.7.  For  the  how 
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Figure  8.4:  Exponent  7  for  the  network  strength  distribution  P{s)  ~  s~'^  plotted  over 
u*{t)l*{t)  with  different  initial  Reynolds  numbers  (green:  Re{to)  =  75,  purple:  Re(to)  =  439, 
yellow:  Reit^)  =  814,  red:  Re{to)  =  1607,  and  blue:  Re{to)  =  2485).  Scale-free  distributions 
are  observed  with  7  coalescing  to  7cr  ~  2.7  up  until  a  bifurcation  at  {u*l*)cr  ~  0.063.  The 
gray  box  shows  7  =  — 2. 7  ±0.5  as  reference. 


under  consideration,  a  bend  in  the  strength  distribution  appears  after  f  30  as  the  system 
starts  to  exhibit  scale  separation.  This  is  caused  by  the  diffusion  of  smaller  scale  structures 
and  their  merging  with  other  structures.  Over  time,  viscous  dissipation  removes  kinetic 
energy  through  the  smaller  eddies  and  leaves  only  the  larger  vortices.  This  behavior  can  be 
described  by  two  power  laws  P{s)  ~  and  P{s)  ~  where  they  capture  the  weaker 
fluid  elements  and  the  larger  stronger  vortices,  respectively.  The  bifurcation  of  these  power 
laws  is  shown  in  Fig.  8.3(c)  indicated  by  the  vertical  dashed  line. 

We  have  considered  a  range  of  Reynolds  numbers  and  observed  that  7  takes  values  of 
7  =  2.7 ±0.5.  The  variations  observed  in  7,  71,  and  72  shown  in  Fig.  8.3(c)  are  influenced  by 
the  chaotic  nature  of  turbulence.  These  parameters  however  appear  to  exhibit  a  coalescing 
behavior  when  they  are  plotted  over  the  product  of  the  characteristic  velocity  and  length, 
u*{t)l*{t).  Here,  we  interpret  u*{t)l*{t)  as  the  circulation  of  vortices  that  have  the  charac¬ 
teristic  velocity  and  length  scales.  As  shown  in  Figure  8.4,  we  observe  that  the  turbulence 
network  shows  coalescence  of  the  scale- free  parameter  7  to  7cr  ~  2.7  over  time  for  different 
cases  of  turbulent  flows.  Once  the  flows  reach  a  state  where  the  characteristic  strength  of 
vortices  is  («*/*)„  ~  0.063,  the  network  distribution  bifurcates  to  display  two  different  slopes 
with  7i  and  72,  as  previously  illustrated  in  Figure  8.3.  This  observation  reveals  that  a  scale- 
free  turbulent  network  is  present  until  the  unforced  turbulent  flow  held  loses  the  smaller-scale 
vortices  and  mostly  contains  vortices  with  strengths  larger  than  the  critical  value  of  («*/*)„• 


8.2  Resilience  of  turbulence  networks 

Characterizing  turbulent  how  with  a  scale-free  network  enables  us  to  view  turbulent  interac¬ 
tions  in  a  systematic  manner  and  provides  insights  into  how  vortical  structures  inhuence  each 
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other.  It  is  known  from  network  analysis  that  scale-free  networks  are  resilient  to  random 
perturbations  but  attacks  towards  network  hubs  can  affect  network  dynamics  in  a  detrimen¬ 
tal  manner  (Albert  et  ai,  2000).  Network  resilience  for  fluid  flow  translates  to  the  difficulty 
of  modifying  the  vortical  interaction  network  and,  consequently,  the  collective  behavior  of 
the  vortices  over  time.  To  measure  the  change  in  vortical  interaction  caused  by  network 
disturbance,  we  can  consider  how  the  removal  of  turbulence  network  nodes  (percolation) 
modihes  the  characteristic  network  length 

^network  =  7  pT  ^  ^min(i(i,j),  (8.3) 

which  is  the  average  shortest  network  distance  d{i,j)  between  any  two  nodes  on  a  network. 
Here,  we  perform  node  percolation  by  setting  the  vorticity  values  at  the  chosen  nodes  to 
be  zero.  The  above  metric  quantifies  how  well  vortical  elements  are  connected  within  a 
turbulent  network.  Note  that  the  distance  here  refers  to  network  distance  based  on  the 
adjacency  matrix  and  not  the  spatial  distance.  In  particular,  we  take  the  inverse  of  each 
adjacency  weight  l/a^  and  evaluate  the  minimal  sum 

+  l/(ikik2  +  •  ■  ■  +  (8-4) 


over  a  network  path  that  connects  nodes  i  and  j  for  this  metric  (Rubinov  &  Sporns,  2010). 
This  metric  /network  can  be  thought  of  as  the  average  of  the  minimal  characteristic  advective 
(commute)  time  per  unit  length  between  every  pair  of  fluid  elements  in  the  domain.  This 
minimal  network  distance  is  determined  using  the  Floyd-Warshall  algorithm  (Floyd,  1962). 

The  changes  in  the  turbulence  network  characteristic  length  /network  when  network  nodes 
are  removed  in  a  random  fashion  and  a  coordinated  manner  targeting  hub  nodes  are  sum¬ 
marized  in  Fig.  8.5.  Here,  the  changes  in  the  normalized  characteristic  network  length 


A/ 


network  — 


^network(^5  ^network(^5  f  0) 

^network(^;  f  0) 


(8.5) 


for  varied  fraction  of  node  removal  are  shown.  While  it  would  be  difficult  to  completely 
remove  nodes  as  we  have  performed  in  this  investigation,  the  present  analysis  sheds  light  on 
how  external  forcing  or  perturbations  can  alter  the  turbulent  flow  from  an  interaction-based 
analysis.  We  observe  that  turbulent  flow  is  resilient  against  random  forcing,  as  evident  from 
the  characteristic  network  length  being  unaffected  even  for  a  large  fraction  /  of  nodes  being 
removed.  This  behavior  is  consistently  observed  over  time.  On  the  other  hand,  we  find  that 
the  global  vortical  interaction  network  can  be  greatly  modified  by  targeting  large  vortex 
cores  (hubs),  as  exhibited  by  the  substantial  change  in  the  characteristic  length.  It  may  be 
more  energetically  expensive  to  remove  well-connected  hub  nodes,  which  often  correspond  to 
regions  of  concentrated  vorticity.  However,  it  is  clear  from  Figure  8.5  that  even  the  smallest 
fraction  of  hub  node  removal  can  influence  the  overall  interaction,  which  suggests  that  hub 
removal  still  provides  a  more  effective  and  efficient  way  to  modify  the  flow  than  random  node 
removal. 
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Figure  8.5:  The  resilience  of  turbulence  network  against  node  removals  for  t  =  15,  30,  and  75 
with  i?e(to)  =  814.  Shown  are  relative  changes  in  the  characteristic  network  length  A/network 
of  turbulent  flow  for  random  node  and  hub  node  removals.  The  colors  of  the  curves  represent 
the  time  when  node  removal  is  considered  and  follows  Figure  8.3. 


When  the  vortical  interaction  network  is  grossly  altered,  the  dynamics  of  the  collection  of 
vortices  would  be  signihcantly  modihed  (Nair  &  Taira,  2015).  These  observations  also  agree 
with  past  studies  in  flow  control  that  identihed  effective  actuation  frequencies  to  be  associated 
with  the  length  scale  of  the  large  coherent  structures  in  turbulent  flows  ( Joslin  &  Miller,  2009; 
Gad-el-Hak,  20006).  With  increasing  time,  we  can  further  notice  that  network  connectivity 
decreases  with  hub  removal  due  to  viscous  dissipation  of  smaller  vortical  structures  and  the 
influence  of  removing  the  core  structures  becomes  more  evident.  The  present  network  based 
understanding  reveals  which  type  of  flow  structures  should  be  targeted  with  flow  control  if 
we  aim  to  alter  the  behavior  of  the  turbulent  flow  held  in  a  global  manner. 
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Chapter  9 

Concluding  Remarks 


Perturbations  of  wall-normal  momentum  and  wall-normal  vorticity  are  injected  into  the 
boundary  layer  to  modify  the  separated  flow  over  a  NACA  0012  airfoil  at  Re  =  23,  000 
and  at  a  high  angle  of  attack.  We  considered  a  massively  separated  baseline  flow  with  a 
separation  region  which  extends  from  the  leading  edge  to  the  trailing  edge.  Eliminating 
the  large  separation  using  actuation  led  to  increased  lift  and  decreased  drag.  In  this  study 
the  control  inputs  were  independently  prescribed  by  a  model  boundary  condition  near  the 
natural  separation  point.  By  examining  the  controlled  flow  field  and  resulting  forces,  we 
quantified  the  necessary  control  inputs  to  the  alter  the  aerodynamic  forces.  In  doing  so,  we 
establish  a  need  for  not  only  a  coefficient  of  momentum  but  also  a  coefficient  of  circulation 
to  account  for  the  angular  momentum  added  to  the  boundary  layer  by  streamwise  vortices. 
The  so-called  coefficient  of  circulation  is  a  function  of  wall-normal  and  azimuthal  velocity, 
accounting  for  the  vorticity  flux  into  boundary  layer. 

Using  the  coefficients  individually  did  not  result  in  an  obvious  trend  that  would  allow 
general  input  arguments  to  quantify  the  resulting  forces.  Therefore,  a  linear  relationship 
is  created  between  the  momentum  directly  added  to  the  flow  (U^)  and  that  added  through 
mixing  the  free  stream  and  boundary  layer  (Ur).  Using  the  total  coefficient  {C^  -|-  Ur) 
identifies  three  regions:  separated  flow,  transition,  and  attached  flow.  Initially,  for  small 
values  of  control  input  lift  decreases  due  to  the  perturbation  added  to  the  boundary  layer, 
but  remains  consistent  for  a  range  of  values  of  Utotai  ^  0.4.  For  these  smaller  inputs  the  size 
of  the  separation  region  remains  similar  to  the  baseline  case.  A  sharp  increase  in  lift  forces 
is  observed  when  the  reverse  flow  region  begins  to  diminish.  Reaching  a  certain  value  for  the 
control  input,  in  this  study  Utotai  ^  0.65%,  flow  is  reattached,  leading  to  improvements  in 
the  aerodynamic  forces.  Realizing  the  total  coefficient  value  necessary  to  reattach  the  flow 
allows  implementation  of  any  actuator  with  Utotai  greater  than  some  critical  value,  which  was 
Utotai  ^  0.65%  for  a  =  9°.  For  example,  we  successfully  identihed  the  value  of  wall-normal 
velocity  (U^j)  that  would  reattach  the  flow  at  ct  =  9°  based  on  the  total  coefficient. 

Moreover,  the  study  has  developed  advanced  analysis  techniques.  First,  the  capability 
to  perform  bi-global  stability  analysis  has  been  developed  and  validated,  which  can  serve 
as  a  basis  for  physics-based  active  flow  control  guided  by  the  knowledge  of  hydrodynamic 
instabilities.  Second,  as  part  of  modeling  complex  unsteady  flows  in  general,  efforts  in 
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this  study  has  led  to  the  initial  development  of  a  novel  approach  in  quantifying  nonlinear 
interactions  present  in  vortical  flows. 

Graph  representations  of  a  system  of  discrete  point  vortices  provide  us  with  a  framework 
for  identifying  important  connections  amongst  vortices.  The  edge  weights  of  the  graph  were 
related  to  the  importance  of  the  ties  between  the  vortices,  which  were  quantihed  based  on  the 
strength  of  the  vortices  and  the  relative  distance  between  them.  We  sparsify  the  connections 
between  these  vortices  based  on  spectral  graph  theory  not  only  to  cut  some  of  these  ties  but 
also  to  redistribute  their  weights  to  preserve  the  dynamics  of  the  vortices.  We  have  demon¬ 
strated  that  an  effective  resistance  approach  to  graph  sparsihcation  proves  extremely  useful 
in  constructing  sparsihed-dynamics  models.  Moreover,  sparsihcation  reduces  the  number 
of  connections  to  0{N  \ogN / e^)  based  on  approximation  order  e.  We  also  observe  an  in¬ 
crease  in  the  effectiveness  of  sparsihcation  with  increasing  number  of  point  vortices.  In  spite 
of  the  nonlinear  nature  of  the  dynamics  of  these  point  vortices,  the  sparse  representations 
are  capable  of  identifying  the  core  vortex  structures  as  well  as  tracking  the  bulk  motion  of 
these  structures  in  an  ehective  manner.  We  have  referred  to  the  dynamical  model  that  uses 
the  sparsihed  vortex  network  as  the  sparsihed-dynamics  model.  In  addition,  the  sparsihed 
models  preserve  spectral  properties  of  the  original  setup.  Resparsihcation  is  utilized  to  en¬ 
hance  the  prediction  of  the  motion  of  the  vortices  over  time.  We  observe  that  sparsihcation 
conserves  the  invariants  of  discrete  vortices  dynamics  such  as  the  Hamiltonian,  linear  and 
angular  impulse,  and  circulation.  The  objectives  of  sparsihed-dynamics  model  are  aligned 
with  those  of  reduced-order  models  but  these  goals  are  achieved  by  sparsifying  the  connec¬ 
tions  rather  than  deriving  a  model  equation  with  state  variable  of  reduced  dimension.  The 
results  of  graph-theoretic  approach  to  discrete  vortex  dynamics  point  to  promising  research 
in  its  applicability  to  a  variety  of  problems  in  huid  mechanics. 

Also,  we  have  identihed  that  the  vortical  interactions  in  two-dimensional  isotropic  tur¬ 
bulence  have  a  scale-free  network  structure  (Taira  et  ai,  2016).  We  have  been  able  to  reveal 
the  structure  by  taking  a  continuous  representation  of  the  how  held  and  quantifying  the  net¬ 
work  using  a  Cartesian  discretization.  For  two-dimensional  isotropic  turbulence,  the  node 
strength  distribution  was  uncovered  to  be  P{s)  ~  s~"',  where  7  =  2.7±0.5.  Furthermore,  we 
have  found  that  the  unforced  turbulent  how  held  possesses  an  underlying  scale-free  network 
structure  until  the  circulation  of  vortices  with  characteristic  velocity  and  length  scales  reach 
{u*l*)cr  ~  0.063.  By  noticing  that  the  turbulence  network  has  scale-free  characteristics,  we 
were  able  to  systematically  show  that  the  turbulence  network  is  resilient  against  random 
perturbations  but  vulnerable  against  coordinated  forcing  on  the  hub  vortices.  It  should  be 
noted  that  estimating  and  controlling  each  and  every  vortical  structure  in  a  turbulent  how  is 
most  likely  improbable  and  impractical.  Instead,  network  analysis  may  provide  a  refreshing 
view  point  on  how  one  can  predict  and  modify  the  collective  dynamics  of  vortices  in  the 
turbulent  how  helds.  We  believe  that  the  network-based  analysis  and  control  (Mesbahi  & 
Egerstedt,  2010;  Liu  et  al,  2011;  Cornelius  et  al,  2013;  Kaiser  et  a/.,  2014&;  Yan  et  al,  2015) 
will  provide  novel  mathematical  fabric  for  paving  the  path  towards  network-based  modeling 
and  control  of  turbulent  hows,  which  can  potentially  impact  a  wide  spectrum  of  problems. 
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Appendix  A 
Control  Inpnt 


A.l  a  =  6° 

The  velocities  used  to  control  the  flow  at  a  =  6°  and  resulting  total  coefficient.  The  accom¬ 
panying  lift  and  drag  coefficients  are  listed  as  well. 


Case 

"^n^max/UoQ 

'^O^raax/UoQ 

C'total[%] 

Rot.  dir. 

Cd 

Cl 

Baseline 

- 

- 

- 

- 

0.061 

0.612 

6A 

0 

1.26 

0.010 

CTR 

0.056 

0.643 

0 

5.05 

0.040 

CTR 

0.043 

0.623 

6B 

0.32 

0 

0.082 

- 

0.047 

0.614 

0.32 

1.26 

0.325 

ROT 

0.043 

0.608 

6C 

0.32 

1.26 

0.325 

CTR 

0.047 

0.615 

0.63 

0 

0.272 

- 

0.050 

0.595 

0.63 

1.26 

0.658 

ROT 

0.043 

0.609 

0.63 

1.26 

0.658 

CTR 

0.049 

0.600 

1.26 

0 

0.911 

- 

0.047 

0.540 

1.26 

1.26 

1.524 

ROT 

0.055 

0.534 

1.26 

1.26 

1.524 

CTR 

0.048 

0.533 

Table  A.l:  All  of  the  control  parameters  used  to  control  for  flow  over  a  NACA0012  airfoil 
at  a  =  6°.  The  specific  cases  used  as  examples  are  denoted. 
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A.2  a  =  9° 


The  velocities  used  to  control  the  flow  at  a  =  9°  and  resulting  total  coefficient.  The  accom¬ 
panying  lift  and  drag  coefficients  are  listed  as  well.  Each  table  contains  different  data.  Table 
A.2  contains  the  majority  of  cases  run,  and  Table  A. 3  contains  the  result  from  modifying 
the  spanwise  spacing  of  the  actuator. 


Case 

'^0  ^max  /  U QQ 

C'total[%] 

Rot.  dir. 

Cd 

Cl 

Baseline 

- 

- 

- 

- 

0.118 

0.584 

0 

5.04 

0.027 

CTR 

0.118 

0.552 

0.631 

0 

0.272 

- 

0.108 

0.416 

0.631 

1.26 

0.657 

CTR 

0.120 

0.411 

0.631 

1.26 

0.657 

ROT 

0.112 

0.414 

0.789 

0.95 

0.736 

CTR 

0.108 

0.429 

0.789 

1.26 

0.848 

CTR 

0.078 

0.679 

0.789 

2.52 

1.294 

CTR 

0.097 

0.609 

0.946 

0 

0.551 

- 

0.118 

0.467 

0.946 

0.946 

0.930 

CTR 

0.101 

0.463 

0.946 

1.26 

1.056 

ROT 

0.065 

0.716 

9E 

0.946 

1.26 

1.056 

CTR 

0.064 

0.726 

0.946 

2.52 

1.561 

CTR 

0.092 

0.628 

9A 

1.26 

0 

0.911 

- 

0.108 

0.403 

1.26 

0.631 

1.217 

ROT 

0.099 

0.493 

1.26 

0.631 

1.217 

CTR 

0.096 

0.510 

9B 

1.26 

0.95 

1.370 

CTR 

0.094 

0.483 

9C 

1.26 

1.26 

1.523 

ROT 

0.069 

0.665 

1.26 

1.26 

1.523 

CTR 

0.067 

0.699 

1.26 

2.523 

2.134 

CTR 

0.069 

0.725 

1.545 

0 

1.747 

- 

0.106 

0.391 

9F 

1.829 

0 

1.298 

— 

0.084 

0.673 

Table  A.2:  Control  inputs  for  flow  over  a  NACA0012  airfoil  at  a  =  9°,  C/c  =  0.1,  and 
i/c  =  0.1.  The  specific  cases  used  as  examples  are  denoted. 
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Iz/  C. 

Case 

^n^max  /  ^oo 

'^O^max/UoQ 

aotal[%] 

Rot.  dir. 

Cd 

Cl 

0.2 

0.892 

0 

0.249 

- 

0.112 

0.457 

0.2 

0.892 

0.892 

0.420 

ROT 

0.113 

0.451 

0.2 

0.892 

1.784 

0.592 

ROT 

0.103 

0.376 

0.2 

1.34 

0 

0.505 

- 

0.114 

0.464 

0.2 

1.34 

0.892 

0.730 

ROT 

0.111 

0.411 

0.2 

1.34 

1.784 

0.955 

ROT 

0.064 

0.731 

0.2 

9F 

1.784 

0 

0.836 

- 

0.075 

0.724 

0.2 

1.784 

1.784 

1.382 

ROT 

0.082 

0.687 

0.067 

1.029 

0 

0.956 

- 

0.100 

0.454 

0.067 

1.029 

1.023 

1.609 

ROT 

0.101 

0.444 

0.05 

0.981 

0 

0.992 

- 

0.096 

0.494 

0.05 

0.891 

0.891 

1.677 

ROT 

0.105 

0.378 

Table  A. 3:  Control  inputs  for  flow  over  a  NACA0012  airfoil  at  a  =  9°,  Izjc  ^  0.1,  and 
i/c  =  0.1.  The  specific  cases  used  as  examples  are  denoted. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Bibliography 


Abe,  Y.,  Okada,  K.,  Sato,  M.,  Nonomura,  T.  &  Fujii,  K.  2013  Significance  of 
three-dimensional  unsteady  flows  inside  of  the  cavity  on  separated-flow  control  around  an 
NACA  0015  using  a  synthetic  jet.  AIAA  2013-2748  . 

Ahuja,  S.  &  Rowley,  C.  W.  2010  Feedback  control  of  unstable  steady  states  of  flow  past 
a  flat  plate  using  reduced-order  estimators.  J.  Fluid  Mech.  645,  447-478. 

Akervik,  E.,  Brandt,  L.,  Henningson,  D.  S.,  Jerome,  H.,  Marxen,  O.  &  Schlat¬ 
ter,  P.  2006  Steady  solutions  of  the  navier-stokes  equations  by  selective  frequency  damp¬ 
ing.  Physics  of  Fluids  18. 

Albert,  R.,  Jeong,  H.  &  Barabasi,  A.-L.  2000  Error  and  attack  tolerance  of  complex 
networks.  Nature  406,  378-382. 

Ashill,  P.  R.,  Flflker,  J.  L.  &  Hackett,  K.  C.  2002  Studies  of  flows  induced  by  sub 
boundary  layer  vortex  generators.  AIAA  Paper  2002-0968  . 

Bagheri,  S.,  Hcepeener,  J.,  Schmid,  P.  &  Henningson,  D.  2009  Input-output  anal¬ 
ysis  and  control  design  applied  to  a  linear  model  of  spatially  developing  flows.  Applied 
Mechanics  Reviews  62  (2). 

Barabasi,  A.-L.  &  Albert,  R.  1999  Emergence  of  scaling  in  random  networks.  Science 
286,  509-512. 

Barrat,  a.,  Barthelemy,  M.  &  Vespignani,  A.  2004  Weighted  evolving  networks: 
coupling  topology  and  weighted  dynamics.  Phys.  Rev.  Let.  92  (22),  228701. 

Barrett,  R.,  Berry,  M.,  Chan,  T.,  Demmel,  J.,  Donato,  J.,  Dongarra,  V.  E., 
Pozo,  R.,  Romaine,  C.  &  vandenHorst,  H.  1993  Templates  for  the  solutions  of 
linear  systems:  Building  blocks  for  iterative  methods.  In  Philidelphia.  SIAM. 

Batchelor,  G.  2000  An  introduction  to  fluid  dynamics.  Cambridge  university  press. 

Benczur,  a.  &  Karger,  D.  R.  1996  Approximating  s  —  t  minimum  cuts  in  0{rR)  time. 
In  Proceedings  of  the  twenty-eighth  annual  ACM  symposium  on  Theory  of  computing,  pp. 
47-55.  ACM. 


77 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Benzi,  R.  &  CoLELLA,  M.  1992  A  simple  point  vortex  model  for  two-dimensional  decaying 
turbulence.  Phys.  Fluids  A  4  (5),  1036-1039. 

Benzi,  R.,  Paladin,  G.  &  Vulpiani,  A.  1990  Power  spectra  in  two-dimensional  turbu¬ 
lence.  Phys.  Rev.  A  42  (6),  3654-3656. 

Berkooz,  G.,  Holmes,  P.  &  Lumley,  J.  1993  The  proper  orthogonal  decomposition  in 
the  analysis  of  turbulent  flows.  Annual  review  of  fluid  mechanics  25  (1),  539-575. 

Boffetta,  G.  &  Ecke,  R.  E.  2012  Two-dimensional  turbulence.  Annu.  Rev.  Fluid  Mech. 
44,  427-451. 

Bollobas,  B.  1998  Modern  graph  theory.  Springer. 

Brehm,  G.  &  Easel,  H.  F.  2011  An  initial  value  problem  approach  to  investigate  biglobal 
stability  problems.  AIAA. 

Brunton,  S.  L.  &  Noack,  B.  R.  2015  Closed-loop  turbulence  control;  progress  and 
challenges.  App.  Mech.  Rev.  67  (5),  050801. 

Galdarelli,  G.  2007  Scale-free  networks.  Oxford  Univ.  Press. 

Ganuto,  G.,  Hussaini,  M.  Y.,  Quarteroni,  A.  &  Zang,  T.  A.  1988  Spectral  methods 
in  fluid  dynamics.  Springer- Verlag,  New  York. 

Gattafesta,  L.,  Song,  Q.,  Williams,  D.,  Rowley,  G.  &  Alvi,  F.  2008  Active  control 
of  flow-induced  cavity  oscillations.  Progress  In  Aerospace  Sciences  44,  479-502. 

Gattafesta,  L.  N.  &  Sheplak,  M.  2011  Actuators  for  active  flow  control.  Annu.  Rev. 
Fluid  Mech.  43,  247-272. 

Gaughemeza,  S.,  a,  B.,  Marghbanks,  T.  L.,  Fagan,  R.  P.,  S.  Ostroff,  N.  M.  F., 
SWERDLOW,  D.  &  WORKING  GROUP,  P.  H.  2011  Role  of  social  networks  in  shaping 
disease  transmission  during  a  community  outbreak  of  2009  HlNl  pandemic  influenza. 
Proceedings  of  the  National  Academy  of  Sciences  108  (7),  2825-2830. 

Ghang,  P.  K.  1976  Control  of  Flow  Separation.  Hemisphere  Publishingn  Corporation. 

Chen,  W.  K.  2004  The  Electrical  Engineering  Handbook.  Elsevier  Science. 

Choi,  H.,  Moin,  P.  &  Kim,  J.  1993  Direct  numerical  simulation  of  turbulent  flow  over 
riblets.  Journal  of  Eluid  Mechanics  255,  503-539. 

Cohen,  R.  &  Havlin,  S.  2010  Complex  Networks:  Structure,  Robustness  and  Eunction. 
Cambridge  Univ.  Press. 

Compton,  D.  A.  &  Johnston,  J.  P.  1992  Streamwise  vortex  production  by  pitched  and 
skewed  jets  in  a  turbulent  boundary  layer.  AIAA  Journal  30,  640-647. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


CoRKE,  T.  C.,  Enloe,  C.  L.  &  Wilkinson,  S.  P.  2010  Dielectric  barrier  discharge 
plasma  actuators  for  flow  control.  Annual  Review  of  Fluid  Mechanics  42,  505-529. 

Cornelius,  S.  P.,  Kath,  W.  K.  &  Motter,  A.  E.  2013  Realistic  control  of  network 
dynamics.  Nature  Comm.  . 

COTTET,  G.-H.  &  Koumoutsakos,  P.  D.  2000  Vortex  methods:  theory  and  practice. 
Cambridge  Univ.  Press. 

Davidson,  P.  A.  2004  Turbulence:  an  introduction  for  scientists  and  engineers.  Oxford 
University  Press. 

Deng,  S.,  Jiang,  L.  &  Liu,  C.  2007  DNS  for  flow  separation  control  around  an  airfoil  by 
pulsed  jets.  Computers  and  Fluids  36  (6),  1040-1060. 

Dorogovtsev,  S.  N.  2010  Lectures  on  Complex  Networks.  Oxford  Univ.  Press. 

Duarte- Cavajalino,  J.  M.,  Jahanshad,  N.,  Lenglet,  C.,  MgMahon,  K.  L., 
DE  ZUBIGARAY,  G.  I.,  MARTIN,  N.  G.,  WRIGHT,  M.  J.,  THOMPSON,  P.  M.  & 
Sapiro,  G.  2012  Hierarchical  topological  network  analysis  of  anatomical  human  brain 
connectivity  and  differences  related  to  sex  and  kinship.  Neuroimage  59,  3784-3804. 

Ellens,  W.,  Spieksma,  F.,  Van  Mieghem,  P.,  Jamakovig,  A.  &  Koou,  R.  2011 
Effective  graph  resistance.  Linear  algebra  and  its  applications  435  (10),  2491-2506. 

Farazmand,  M.  M.,  Kevlahan,  N.  K.-R.  &  Protas,  B.  2011  Controlling  the  dual 
cascade  of  two-dimensional  turbulence.  J.  Fluid  Mech.  668,  202-222. 

Floyd,  R.  W.  1962  Algorithm  97:  Shortest  path.  Comm.  ACM  5  (6),  345. 

Frisgh,  U.  1995  Turbulence.  Cambridge  Univ.  Press. 

Gad-el-Hak,  M.  2000a  Flow  control:  passive,  active,  and  reactive  flow  management.  Cam¬ 
bridge  Univ.  Press,  London. 

Gad-el-Hak,  M.  2000b  Flow  control:  passive,  active,  and  reactive  flow  management.  Cam¬ 
bridge  Univ.  Press. 

Gilarranz,  j.  L.,  Traub,  L.  W.  &  Redinoitis,  O.  K.  2005  A  new  class  of  synthetic 
jet  actuators-Part  II:  application  to  flow  separation  control.  Journal  of  Fluids  Engineering 

127. 

Glass,  R.  J.,  Glass,  L.  M.,  Beyeler,  W.  E.  &  Min,  H.  J.  2006  Targeted  social 
distancing  design  for  pandemic  influenza.  Emerg.  Infect.  Diseases  12  (11),  1671-1681. 

Glezer,  a.  &  Amitay,  M.  2002  Synthetic  jets.  Annu.  Rev.  Fluid  Mech.  34,  503-529. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Greenblatt,  D.,  Wygnanski,  I.  J.  &  Rumsey,  C.  L.  2015  Aerodynamic  flow  control. 
Encyclopedia  of  Aerospace  Engineering  . 

Greengard,  L.  &  Rokhlin,  V.  1987  A  fast  algorithm  for  particle  summations.  J.  Comput. 
Phys.  73,  325-348. 

Ham,  F.  &  Iaggarino,  G.  2004  Energy  conservation  in  collocated  discretization  schemes 
on  unstructured  meshes.  Annual  research  brief.  Center  for  Turbulence  Research,  Stanford 
University. 

Ham,  F.,  Mattson,  K.  &  Iaggarino,  G.  2006  Accurate  and  stable  finite  volume  opera¬ 
tors  for  unstructured  flow  solvers.  Tech.  Rep..  Center  for  Turbulence  Research. 

Hemati,  M.,  Eldredge,  J.  D.  &  Speyer,  J.  L.  2014  Improving  vortex  models  via 
optimal  control  theory.  Journal  of  Eluids  and  Structures  . 

Hinze,  j.  O.  1975  Turbulence.  McGraw-Hill,  New  York. 

Holmes,  P.,  Lumley,  J.  L.  &  Berkooz,  G.  1996  Turbulence,  coherent  structures,  dy¬ 
namical  systems  and  symmetry.  Cambridge  Univ.  Press. 

Hornung,  H.  1989  Vorticity  generation  and  transport.  In  Tenth  Australasian  Eluid  Me¬ 
chanics  Conference  -  University  of  Melbourne. 

Huang,  L.,  Huang,  P.  G.  &  LeBeau,  R.  P.  2004  Numerical  study  of  blowing  and 
suction  control  mechanism  on  NACA  0012  airfoil.  Journal  of  Aircraft  41  (1). 

JosLiN,  R.  D.  &  Miller,  D.  (ed.)  2009  Eundamentals  and  applications  of  modern  flow 
control.  AIAA. 

Kaiser,  E.,  Noagk,  B.  R.,  Gordier,  L.,  Spohn,  A.,  Segond,  M.,  Abel,  M., 

Daviller,  G.  &  Niven,  R.  K.  2014a  Cluster-based  reduced-order  modelling  of  a  mixing 
layer.  Journal  of  Eluid  Mechanics  754,  365-414. 

Kaiser,  E.,  Noagk,  B.  R.,  Gordier,  L.,  Spohn,  A.,  Segond,  M.,  Abel,  M., 

Daviller,  G.,  Osth,  J.,  Krajnovic,  S.  &  Niven,  R.  K.  2014&  Cluster-based 

reduced-order  modelling  of  a  mixing  layer.  J.  Eluid  Mech.  754,  365-414. 

Keener,  J.  &  Levin,  A.  2011  Spectral  sparsiflcation  in  the  semi-streaming  setting.  Leibniz 
International  Proceedings  in  Informatics  (LIPIcs)  series  9,  440-451. 

Kerho,  M.,  Hutgherson,  S.,  Blagkwelder,  R.  F.  &  H.,  L.  R.  1993  Vortex  genera¬ 
tors  used  to  control  laminar  separation  bubbles.  Journal  of  Aircraft  30  (3),  315-319. 

Kitsios,  V.,  Rodriguez,  D.,  Theofilis,  V.,  Ooi,  A.  &  Soria,  J.  2009  Biglobal 
stability  analysis  in  curvilinear  coordinates  of  massively  separted  lifting  bodies.  Journal 
of  Computational  Physics  228,  7181-7196. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Klein,  D.  J.  &  Randic,  M.  1993  Resistance  distance.  Journal  of  Mathematical  Chemistry 
12  (1),  81-95. 

Kojima,  R.,  Nonomura,  T.,  Oyama,  A.  &  Fujii,  K.  2013  Large-eddy  simnlation  of 
low-Reynolds-number  flow  over  thick  and  thin  NACA  airfoils.  Journal  of  Aircraft  50  (1), 
187-196. 

Kotapati,  R.  B.,  Mittal,  R.,  Marxen,  O.,  Ham,  F.,  You,  D.  &  Cattafesta,  L.  N. 
2010  Nonlinear  dynamics  and  synthetic-jet-based  control  of  a  canonical  separated  flow.  J. 
Fluid  Meek  654,  65-97. 

Kraichnan,  R.  H.  &  Montgomery,  D.  1980  Two-dimensional  tnrbulence.  Rep.  Prog. 
Rhys.  43,  547-619. 

Lachmann,  G.  V.  1961  Boundary  layer  and  flow  control.  Its  principles  and  applications, 
vol  1  &  2.  Pergamon  Press. 

Lehoucq,  R.  B.,  Sorensen,  D.  C.  &  Yang,  C.  1998  Arpack  nsers  guide:  Solution  of 
large-scale  eigenvalue  problems  with  implicitly  restarted  arnoldi  methods.  SIAM. 

Leonard,  A.  1980  Vortex  methods  for  flow  simulation.  Journal  of  Computational  Physics 
37  (3),  289-335. 

Lesieur,  M.  2008  Turbulence  in  fluids,  4th  edn.  Springer. 

Lin,  j.  C.  2002  Review  of  research  on  low-prohle  vortex  generators  to  control  boundary-layer 
separation.  Progress  in  Aerospace  Science  38,  389-420. 

Lin,  j.  C.,  Robinson,  S.  K.,  Mgghee,  R.  J.  &  Valarezo,  W.  O.  1994  Separation 
control  on  hihg-lift  airfoil  via  micro- vortex  generators.  Journal  of  Aircraft  31  (6). 

Little,  J.,  Nishihara,  M.,  Adamovigh,  I.  &  Samimy,  M.  2010  High-lift  airfoil  trail¬ 
ing  edge  separation  control  using  a  single  dielectric  barrier  discharge  plasma  actuator. 
Experiments  in  Fluids  48,  521-537. 

Liu,  Y.-Y.,  Slotine,  J.-J.  &  Barabasi,  A.-L.  2011  Controllability  of  complex  networks. 
Nature  473  (7346),  167-173. 

Lloyd-Smith,  j.,  Sghreiber,  S.,  Kopp,  P.  &  Getz,  W.  2005  Superspreading  and  the 
effect  of  individual  variation  on  disease  emergence.  Nature  438,  355-359. 

Ma,  Z.,  Ahuja,  S.  &  Rowley,  C.  W.  2009  Reduced  order  models  for  control  of  fluids 
using  the  eigensystem  realization  algorithm.  Theo.  Comp.  Fluid  Dyn.  25  (1),  233-247. 

Mg  Williams,  J.  C.  1984  The  emergence  of  isolated  coherent  vortices  in  turbulent  flow.  J. 
Fluid  Mech.  146,  21-43. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Mesbahi,  M.  &  Egerstedt,  M.  2010  Graph  theoretic  methods  in  multiagent  networks. 
Princeton  Univ.  Press. 

Mieghem,  P.  V.  2011  Graph  spectra  for  complex  networks.  Cambridge  University  Press. 

Mohar,  B.  1991  The  Laplacian  spectrum  of  graphs.  In  Graph  theory,  combinatorics,  and 
applications  (ed.  Y.  Alavi,  G.  Chartrand,  O.  Ollermann  &  A.  Schwenk),  pp.  871-898. 
Wiley,  New  York. 

Morinishi,  Y.,  Lund,  T.  S.,  Vasilyev,  O.  V.  &  Mom,  P.  1998  Fully  conservative  high 
order  hnite  difference  schemes  for  incompressible  flow.  Journal  of  Gomputational  Physics 
143,  90-124. 

Morris,  M.  1993  Epidemiology  and  social  networks  -  modeling  structured  diffusion.  Sociol. 
Method  Res.  22,  99-126. 

Munday,  P.  M.  &  Taira,  K.  2014  Separation  control  on  NACA  0012  airfoil  using  mo¬ 
mentum  and  wall-normal  vorticity  injection.  AIAA  2014-2685  . 

Munday,  P.  M.  &  Taira,  K.  2015  Surface  vorticity  flux  analysis  in  separation  control  on 
NACA  0012  airfoil.  AIAA  2015-2632  . 

Muppidi,  S.  &  Mahesh,  K.  2005  Study  of  trajectories  of  jets  in  crossflow  using  direct 
numerical  simulations.  Journal  of  Fluid  Mechanics  530,  81-100. 

Nair,  a.  G.  &  Taira,  K.  2015  Network-theoretic  approach  to  sparsihed  discrete  vortex 
dynamics.  J.  Fluid  Mech.  768,  549-571. 

Natarajan,  R.  &  Agrivos,  A.  1993  The  instability  of  the  steady  flow  past  spheres  and 
disks.  Journal  of  Fluid  Mechanics  254,  323-344. 

Newman,  M.  E.  J.  2004  Fast  algorithm  for  detecting  community  structure  in  networks. 
Physical  review  E  69,  066133. 

Newman,  M.  E.  J.  2010  Networks:  an  introduction.  Oxford  Univ.  Press. 

Newton,  P.  K.  2001  The  N -vortex  problem:  analytical  techniques,  Applied  Mathematical 
Sciences,  vol.  145.  Springer. 

Noagk,  B.  R.  &  Egkelmann,  H.  1994  A  global  stability  analysis  of  the  steady  and 
periodic  cylinder  wake.  Journal  of  Fluid  Mechanics  270,  297-330. 

Noagk,  B.  R.,  Papas,  P.  &  Monkewitz,  P.  A.  2005  The  need  for  a  pressure-term 
representation  in  empirical  Galerkin  models  of  incompressible  shear  flows.  J.  Fluid  Mech. 
523,  339-365. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Owen,  J.  P.,  Li,  Y.-O.,  Ziv,  E.,  Strominger,  Z.,  Gold,  J.,  Bukhpun,  P., 
Wakahiro,  M.,  Friedman,  E.  J.,  Sherr,  E.  H.  &  Mukherjee,  P.  2013  The  struc¬ 
tural  connectome  of  the  human  brain  in  agenesis  of  the  corpus  callosum.  Neuroimage  70, 
340-355. 

Pedro,  H.  T.  C.  &  Kobayashi,  M.  H.  2008  Numerical  study  of  stall  delay  on  humpback 
whale  flippers.  AIAA  Paper  2008-0584. 

Peleg,  D.  &  Ullman,  j.  1989  An  optimal  synchronizer  for  the  hypercube.  SIAM  Journal 
on  computing  18  (4),  740-747. 

Pope,  S.  B.  2000  Turbulent  flows .  Cambridge  Univ.  Press. 

Porter,  M.  A.,  Mugha,  P.  J.,  Newman,  M.  E.  J.  &  Warmbrand,  C.  M.  2005  A 
network  analysis  of  committees  in  the  U.S.  House  of  Representatives.  Proc.  Nat.  Acad. 
Sci.  102  (20),  7057-7062. 

Raju,  R.,  Aram,  E.,  Mittal,  R.  &  Cattafesta,  L.  2009  Simple  models  of  zero-net 
mass-flux  jets  for  flow  control  simulations.  International  Journal  of  Flow  Control  1  (3), 
179-197. 

Rathay,  N.,  Bougher,  M.,  Amitay,  M.  &  Whalen,  E.  2014a  Parametric  study  of 
synthetic-jet-based  control  for  performance  enhancement  of  a  vertical  tail.  AIAA  Journal 
52  (11),  2440-2454. 

Rathay,  N.,  Bougher,  M.,  Amitay,  M.  &  Whalen,  E.  20146  Performance  enhance¬ 
ment  of  a  vertical  tail  using  synthetic  jet  actuators.  AIAA  Journal  52  (4). 

Robinson,  K.,  Cohen,  T.  &  Colijn,  C.  2012  The  dynamics  of  sexual  contact  networks: 
effects  on  disease  spread  and  control.  Theo.  Popul.  Bio.  81,  89-96. 

Rodriguez,  D.  2010  Global  Instability  of  Laminar  Separation  Bubbles.  PhD  Thesis,  Uni- 
versidad  Politecnica  De  Madrid. 

Rowley,  C.,  Colonies,  T.  &  Murray,  R.  2004  Model  reduction  for  compressible  flows 
using  POD  and  Galerkin  projection.  Physica  D:  Nonlinear  Phenomena  189  (1),  115-129. 

Rowley,  C.  W.  2005  Model  reduction  for  fluids,  using  balanced  proper  orthogonal  decom¬ 
position.  Int.  J.  Bif.  Chaos  15  (3),  997-1013. 

Rowley,  C.  W.,  Mezic,  I.,  Bagheri,  S.  &  Henningson,  D.  S.  2009  Spectral  analysis 
of  nonlinear  flows.  J.  Fluid  Mech.  641,  115-127. 

Rubinov,  M.  &  Sporns,  O.  2010  Complex  network  measures  of  brain  connectivity:  Uses 
and  interpretations.  Neuroimage  52,  1059-1069. 

Saffman,  P.  G.  1992  Vortex  dynamics.  Cambridge  Univ.  Press. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Salathe,  M.  &  Jones,  J.  H.  2010  Dynamics  and  control  of  diseases  in  networks  with 
commnnity  strncture.  PLoS  Comput.  Bio.  6  (4),  el000736. 

Sato,  M.,  Aono,  H.,  Yakeno,  A.,  Nonomura,  T.,  Fuji:,  K.,  Okada,  K.  &  Asada, 
K.  2015  Mnltifactorial  effects  of  operating  conditions  of  dielectric-barrier-discharge  plasma 
actnator  on  laminar-separated-flow  control.  AIAA  Journal  53  (9). 

Schmid,  P.  J.  2007  Nonmodal  stability  theory.  Annu.  Rev.  Fluid  Mech.  39,  129-162. 

Schmid,  P.  J.  2010  Dynamic  mode  decomposition  of  nnmerical  and  experimental  data. 
Journal  of  Fluid  Mechanics  656,  5-28. 

Seele,  R.,  Graff,  E.,  Gharib,  M.,  Taubert,  L.,  Lin,  J.  &  Wygnanski,  I.  2012 
Improving  rndder  effectiveness  with  sweeping  jet  actnators.  AIAA  paper  3244. 

Seele,  R.,  Graff,  E.,  Lin,  J.  &  Wygnanski,  I.  2013  Performance  enhancement  of  a 
vertical  tail  model  with  sweeping  jet  actnators.  AIAA  Paper  411,  2013. 

Seifert,  A.  &  Pack,  L.  T.  1999  Oscillatory  excitation  of  nnsteady  compressible  flows 
over  airfoils  at  flight  reynolds  nnmbers.  AIAA  10.2514/6.1999-925  . 

Selby,  G.  V.,  Lin,  J.  G.  &  Howard,  F.  G.  1992  Control  of  low-speed  tnrbulent  separated 
flow  using  jet  vortex  generators.  Experiments  in  Fluids  12,  394-400. 

Shan,  H.,  Jiang,  L.,  Ghaoqun,  L.,  Love,  M.  &  Maines,  B.  2008  Numerical  study 
of  passive  and  active  flow  separation  control  over  a  NACA  0012  airfoil.  Computers  and 
Fluids  pp.  975-992. 

Skillen,  a.,  Revell,  A.,  Pinelli,  A.,  Piomelli,  U.  &  Favier,  J.  2015  Flow  over  a 
wing  with  leading-edge  undulations.  AIAA  Journal  53  (2),  464-472. 

Spielman,  D.  a.  &  Srivastava,  N.  2011  Graph  sparsihcation  by  effective  resistances. 
SIAM  J.  Comput.  40  (6),  1913-1926. 

Spielman,  D.  A.  &  Teng,  S.-H.  2011  Spectral  sparsihcation  of  graphs.  SIAM  J.  Comput. 
40  (4),  981-1025. 

Srivastava,  N.  2010  Spectral  sparsification  and  restricted  invertibility .  PhD  Thesis,  Yale 
University. 

Taira,  K.,  Nair,  A.  G.  &  Brunton,  S.  L.  2016  Network  structure  of  two-dimensional 
decaying  isotropic  turbulence.  J.  Fluid  Mech.  795,  R2. 

Taylor,  G.  I.  1918  On  the  dissipation  of  eddies.  Reports  and  Memoranda  598.  Aero.  Res. 
Comm. 

Tennekes,  H.  &  Lumley,  J.  L.  1972  A  first  course  in  turbulence.  MIT  Press,  Cambridge. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Theofilis,  V.  2003  Advances  in  global  linear  instability  analysis  of  nonparallel  and  three- 
dimensional  flows.  Progress  In  Aerospace  Sciences  39,  249-315. 

Theofilis,  V.  2011  Global  linear  instability.  Annual  Review  of  Fluid  Mechanics  43,  319- 
352. 

Theofilis,  V.,  Duck,  P.  W.  &  Owen,  J.  2004  Viscous  linear  stability  analysis  of  rect¬ 
angular  duct  and  cavity  flows.  Journal  of  Fluid  Mechanics  505,  249-286. 

Vreman,  a.  W.  2004  An  eddy- viscocity  subgrid-scale  model  for  turbulent  shear  flow: 
Algebraic  theory  and  applications.  Physics  of  Fluids  16  (10). 

Wang,  C.  &  Eldredge,  J.  2013  Low-order  phenomenological  modeling  of  leading-edge 
vortex  formation.  Theoretical  and  Computational  Fluid  Dynamics  27  (5),  577-598. 

Weiss,  J.  B.,  Provenzale,  A.  &  MgWilliams,  J.  C.  1998  Lagrangian  dynamics  in 
high-dimensional  point-vortex  systems.  Phys.  Fluids  10  (8),  1929-1941. 

Whalen,  E.  A.,  Lagy,  D.  S.,  Lin,  J.  C.,  Andino,  M.  Y.,  Washburn,  A.  E.,  Graff, 
E.  C.  &  Wygnanski,  I.  J.  2015  Performance  enhancement  of  a  full-scale  vertical  tail 
model  equipped  with  active  flow  control.  In  53rd  AIAA  Aerospace  Sciences  Meeting,  pp. 
1-11. 

Wu,  J.-Z.,  Lu,  X.-Y.,  Denny,  A.  D.,  Fan,  M.  &  Wu,  J.-M.  1998  Post-stall  flow 
control  on  an  airfoil  by  local  unsteady  forcing.  J.  Fluid  Mech.  371,  21-58. 

Wu,  J.-Z.,  Ma,  H.-Y.  &  Zhou,  M.-D.  2006  Vorticity  and  vortex  dynamics.  Springer- 
Verlag. 

Yan,  G.,  Tsekenis,  G.,  Barzel,  B.,  Slotine,  J.-J.,  Liu,  Y.-Y.  &  Barabasi,  A.-L. 
2015  Spectrum  of  controlling  and  observing  complex  networks.  Nature  Physics  . 

You,  D.,  Ham,  F.  &  Moin,  P.  2008  Discrete  conservation  principles  in  large-eddy  simu¬ 
lation  with  application  to  separation  control  over  an  airfoil.  Physics  of  Fluids  20. 

Zhang,  W.  &  Samtaney,  R.  2016  Biglobal  linear  stability  analysis  on  low-re  flow  past 
an  airfoil  at  high  angle  of  attack.  Physics  of  Fluids  (1994-present)  28  (4),  044105. 

Zhang,  X.  2003  The  evolution  of  co-rotating  vortices  in  a  canonical  boundary  laryer  with 
inclined  jets.  Physics  of  Fluids  15  (12),  3693-3702. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


AhOiSH  Deliverables  Submission  Survey 


Response  ID:6661  Data 


1. 

1 .  Report  Type 
Final  Report 
Primary  Contact  E-mail 

Contact  email  if  there  is  a  problem  with  the  report. 

ktaira@fsu.edu 

Primary  Contact  Phone  Number 

Contact  phone  number  if  there  is  a  problem  with  the  report 

850-645-0140 

Organization  /  Institution  name 

Florida  State  University 

Grant/Contract  Title 

The  full  title  of  the  funded  effort. 

(YIP  1 3)  Understanding  the  fundamental  roles  of  momentum  and  vorticity  injeotions  in  flow  control 

Grant/Contract  Number 

AFOSR  assigned  control  number.  It  must  begin  with  "FA9550"  or  "F49620"  or  "FA2386". 

FA9550-13-1-0183 
Principal  Investigator  Name 

The  full  name  of  the  principal  investigator  on  the  grant  or  contract. 

Kunihiko  Taira 
Program  Manager 

The  AFOSR  Program  Manager  currently  assigned  to  the  award 
Douglas  Smith 
Reporting  Period  Start  Date 

05/1 5/2013 

Reporting  Period  End  Date 
05/1 4/2016 
Abstract 


The  objective  of  this  study  is  to  numerically  investigate  the  fundamental  roles  that  momentum  and  vorticity 
injections  play  in  suppressing  flow  separation  over  a  canonical  airfoil.  Open-loop  control  of  separated, 
incompressible  flow  over  a  NACA  001 2  airfoil  at  Re  =  23,000  is  examined  through  large-eddy  simulations. 
We  find  that  the  modification  to  the  flow  field  can  be  captured  by  quantifying  both  the  effects  of  wall-normal 
momentum  (coefficient  of  momentum)  and  wall-normal  vorticity  (derived  coefficient  of  circulation),  by 
considering  a  newly  defined  total  input  parameter  (total  coefficient).  Moreover,  the  study  has  developed 
advanced  analysis  techniques.  First,  the  capability  to  perform  bi-global  stability  analysis  has  been 
developed  and  validated,  which  can  serve  as  a  basis  for  physics-based  active  flow  control  guided  by  the 
knowledge  of  hydrodynamic  instabilities.  Second,  as  part  of  modeling  complex  unsteady  flows  in  general, 
efforts  in  this  study  have  led  to  the  initial  development  of  a  novel  network-theoretic  approach  in  quantifying 
nonlinear  interactions  present  in  vortical  flows. 


Distribution  Statement 

This  is  biock  12  on  the  SF298  form. 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Distribution  A  -  Approved  for  Public  Release 

Explanation  for  Distribution  Statement 

If  this  is  not  approved  for  pubiic  reiease,  piease  provide  a  short  expianation.  E.g.,  contains  proprietary  information. 
SF298  Form 

Please  attach  your  SF298  form.  A  blank  SF298  can  be  found  here.  Please  do  not  password  protect  or  secure  the  PDF 
The  maximum  fiie  size  for  an  SF298  is  50MB. 

AFOSR-YIP-SF298.pdf 

Upload  the  Report  Document.  File  must  be  a  PDF.  Please  do  not  password  protect  or  secure  the  PDF .  The 
maximum  file  size  for  the  Report  Document  is  50MB. 

AFOSR_YIP_final_report.pdf 

Upload  a  Report  Document,  if  any.  The  maximum  file  size  for  the  Report  Document  is  50MB. 

Archival  Publications  (published)  during  reporting  period: 

P.  M.  Monday  &  K.  Taira,  "Separation  control  on  NACA  001 2  airfoil  using  momentum  and  wall-normal 
vorticity  injection,"  AIAA  2014-2685,  2014. 

P.  M.  Monday  &  K.  Taira,  "Surface  vorticity  flux  analysis  in  separation  control  on  NACA  001 2  airfoil,"  AIAA 
2015-2632,2015. 

A.  G.  Nair  &  K.  Taira,  "Network-theoretic  approach  to  sparsified  discrete  vortex  dynamics,"  J.  Fluid  Mech. 
768,549-571,2015. 

K.  Taira,  A.  G.  Nair,  &  S.  L.  Brunton,  "Network  structure  of  two-dimensional  decaying  isotropic  turbulence," 

J.  Fluid  Mech.  795,  R2,2016. 

2.  New  discoveries,  inventions,  or  patent  disclosures: 

Do  you  have  any  discoveries,  inventions,  or  patent  disclosures  to  report  for  this  period? 

Yes 

Please  describe  and  include  any  notable  dates 

Pending  Patent; 

K.  Taira,  P.  Monday,  and  F.  AIvi,  No.  61/947,164,  Swirling  jet  actuator  for  control  of  separated  and  mixing 
flows,  201 5. 

Do  you  plan  to  pursue  a  claim  for  personal  or  organizational  intellectual  property? 

Yes 

Changes  in  research  objectives  (if  any): 

Change  in  AFOSR  Program  Manager,  if  any: 

Extensions  granted  or  milestones  slipped,  if  any: 

AFOSR  LRIR  Number 
LRIR  Title 
Reporting  Period 
Laboratory  Task  Manager 
Program  Officer 
Research  Objectives 
Technical  Summary 

Funding  Summary  by  Cost  Category  (by  FY,  $K) 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Starting  FY 

FY+1 

FY+2 

Salary 

Equipment/Facilities 

Supplies 

Total 

Report  Document 

Report  Document  -  Text  Analysis 

Report  Document  -  Text  Analysis 

Appendix  Documents 

2.  Thank  You 

E-mail  user 

Aug  1 2,  201 6  1 0:48:54  Success:  Email  Sent  to:  ktaira@fsu.edu 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


